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In this work, we investigate the influence of the anisotropic heating on the spontaneous instability and 
evolution of thin Harris-type collisionless current sheets, embedded in antiparallel magnetic fields. In 
particular, we explore the influence of the macroparticle shape-function using a 2D version of the PIC 
code ACRONYM. We also investigate the role of the numerical collisionality due to the finite number of 
macroparticles in PIC codes. It is shown that it is appropriate to choose higher order shape functions of the 
macroparticles compared to a larger number of macroparticles per cell. This allows to estimate better the 
anisotropic electron heating due to the collisions of macroparticles in a PIC code. Temperature anisotropies 
can stabilize the tearing mode instability and trigger additional current sheet instabilities. We found a 
good agreement between the analytically derived threshold for the stabilization of the anisotropic tearing 
mode and other instabilities, either spontaneously developing or initially triggered ones. Numerical effects 
causing anisotropic heating at electron time scales, become especially important for higher mass ratios 
(above rrii/me = 180). If numerical effects are carefully taken into account, one can recover the theoretical 
estimated linear growth rates of the tearing instability of thin isotropic collisionless current sheets, also for 
higher mass ratios. 
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I. INTRODUCTION 

Current sheets (CS) are ubiquitous in space plasmas, 
as those of the solar corona or of the Earth’s magnetotail. 
They can become unstable and enable the dissipation 
of magnetic energy by reconnection, causing turbulence, 
plasma and particle acceleration^. In fact, CS are the 
preferred locations where magnetic reconnection can take 
place. CS are prone to a number of different instabilities 
driven by currents gradients or streams as source of free 
energy^’^. We focus on thin CS, that have a halfwidth of 
the order of the ion inertial length. This is supposed to be 
the natural limit of their thinning by compression. Inside 
of the ion inertial length, particles may experience res¬ 
onant wave-particle interactions which may provide the 
physical mechanism for dissipation in collisionless plas¬ 
mas"''’^. In fact, many space plasmas in the solar system 
and beyond are collisionless, hence thin CS are of fun¬ 
damental importance in those environments, e.g., for the 
onset of solar flares and substorms®. 

The most fundamental instability of CS is the tear¬ 
ing instability which produces magnetic islands, i.e.: a 
change of topology necessary for magnetic reconnection. 
It was first investigated for collision-dominated plasma 
in the framework of Magnetohydrodynamics (MHD) 
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approach^. Later, it was proposed and investigated also 
for collisionless plasmas®^^^. The mechanism for the col¬ 
lisionless tearing instability and for its non-linear satura¬ 
tion depends on the specific parameters of the CS. It can 
be the lack of gyrotropy of the electron distribution func¬ 
tion, the electron pressure anisotropy and/or the electron 
trapping (see Ref. 12 for a recent review). 

Previous studies have shown that the collisionless tear¬ 
ing instability is very sensitive to electron tempera¬ 
ture anisotropyOn the other hand, it is known 
that numerical simulations by Particle-in-Cell (PIC)^^ 
and Vlasov codes^® may produce artificial heating under 
some circumstances (see Refs. 19 and 20 and references 
therein). Numerical heating may mimic physical heating 
processes and, like those, it may affect the growth rate of 
the tearing instability. In addition, it can cause instabil¬ 
ities that require pre-heating for their development. Two 
independent numerical heating mechanisms were identi¬ 
fied (see Ref. 21 and references therein): grid heating and 
scattering. According to Ref. 21, grid heating is due to 
a kinetic instability resulting from the aliasing of high 
frequency modes not resolved by the grid. Hence, grid 
heating is important when the Debye length is not re¬ 
solved well enough by the grid. Since one can avoid it by 
resolving the Debye length, we will not analyze or discuss 
this well-known process here. Instead, we focus on the 
numerical enhanced particle heating due to scattering. 
This unphysical process arises due to a finite number 
of macroparticles per cell replacing, as an approxima¬ 
tion, the continuous phase space. When a macroparticle 
moves to the next cell, it produces locally random flue- 



2 


tuations of the electric potential that will act backwards 
on the other macroparticles. Thus, macroparticles expe¬ 
rience an effective scattering due to these non-physical 
forces which are absent in real collisionless plasmas (see, 
e.g.. Refs. 22 and 23). Therefore, this scattering can be 
interpreted as a form of numerical collisions among the 
macroparticles. 


II. SIMULATION SETUP 

We consider the stability of CS in a Harris equilib¬ 
rium^^, with a magnetic field changing as 

B{x) = Baoy tanh {x/L) y, (1) 


While in the past extensive studies of the grid scale- 
induced numerical heating and their effects in PIC simu¬ 
lations were carried out, the heating due to numerically 
caused macroparticle scattering remained less well an¬ 
alyzed. Just recently, Cormier-Michel et al.^^ showed 
that there are some effects in laser wakefield accelerators 
depending critically on some threshold, above which nu¬ 
merical scattering dominates those phenomena. Further¬ 
more, these authors found that the choice of an appropri¬ 
ate shape function of macroparticles can be much more 
efficient in reducing the numerical scattering and heat¬ 
ing than the usual approach of increasing the number of 
macroparticles per cell. See Appendix A for the definition 
and importance of the shape function (a form of interpo¬ 
lation scheme) in PIC codes. For a discussion about the 
dependence of the numerical scattering/collisions on the 
shape function, see Appendix B. 


where L is the halfwidth of the CS and B^oy is the asymp¬ 
totic magnetic field. The Harris equilibrium is sustained 
by the current carried by counterstreaming electrons and 
ions (denoted with the subindices e and i) with drifting 
speeds Uoi and Uoe^ respectively. This current density 
is given by: 

J{x) = 2enoUDiCOsh~^ (x/L) z and 

. . (. 2 ) 

The latter condition guarantees the absence of initial 
electric fields. The asymptotic magnetic field and the 
peak density at the center of the CS (with uq = uoe = 
noi) are related by the equilibrium condition between the 
magnetic and thermal pressures 



noiksTe + ksTi). 


( 3 ) 


The relation between shape functions and heating can 
be understood in the following way: the stochastic force 
producing the numerical scattering is introduced by er¬ 
rors in the interpolation of the electromagnetic field at 
the macroparticles position. Thus, it leads to errors in 
the resulting momenta of the macroparticles, such as 
their temperature^^. Therefore, the numerical heating 
will be reduced by smoother interpolations, since they 
will reduce the non-physical forces by smoothing out 
the sharp border of the macroparticles (see Sec. 13.5 of 
Ref. 19). 

So far, the influence of the use of higher order shape 
functions in PIC simulations on the stability of collision¬ 
less Harris CS has not been explicitly investigated yet. 
For this reason, the main purpose of this paper is to study 
the influence of the shape functions on heating processes 
in the course of the CS evolution, as well as to distinguish 
their importance compared to other numerical processes. 
As we will show, if those numerical considerations are 
not taken into account, the simulations may lead to non¬ 
physical results. 

This paper is organized as follows. In Sec. H, we 
describe the simulation setup. In Sec. HI, we discuss 
the consequences of the numerically-induced anisotropic 
heating. In Sec. IV, we analyze the origins and conse¬ 
quences of an initially imposed temperature anisotropy. 
In Sec. V, we analyze the influence of these numerical 
effects on the development of the tearing instability in 
a Harris CS. Finally, we summarize our findings in the 
conclusion. Sec. VI. 


We do not impose a background density in the first place, 
in order to minimize additional effects due to that pop¬ 
ulation and in that way to reduce the number of free 
parameters to be chosen. Note that we will discuss the 
consequences of a background plasma in Sec. HI D 

On purpose, we do not impose an initial perturbation 
other than the shot noise due to the finite number of 
macroparticles, because we focus on the onset and growth 
of the spontaneous instabilities^®’^®. Hence, the system 
evolves only from this macroparticle noise. Note that 
simulations with an initial perturbation bypass this stage 
allowing the system to reach directly the phase of fully 
developed fast reconnection^^. 

See Fig. 1 for a schematics of the setup. 



Figure 1. Schematics of the geometry of the simulation setup. 
The magnetic field B is in the y direction on dependence on 
the X coordinate. Correspondingly, the density n and current 
density J vary along the x direction. The 2D simulations 
are carried out in the x — y plane. The CS is sustained by 
counterstreaming electrons and ions that produce a current 
density J in the out-of-plane direction 2 . 
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The four parameters necessary to fully specify this Har¬ 
ris sheet are chosen according to Ref. 28: 
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( 4 ) 


where c is the speed of light, Wpe = j{e^me) is the 

electron plasma frequency (calculated with the central 
density no), Hce = eBaoy/'me is the electron cyclotron 
frequency in the asymptotic magnetic field. The cho¬ 
sen mass ratio is a compromise between computational 
speed and the possibility to separate the effects of ion 
and electron motion. The width of the CS is of the or¬ 
der of the ion inertial length which is supposed to be 
the natural limit of the thinning of CS. The ratio of fre¬ 
quencies LOpej^ce is choscn to be relatively small to save 
computation time. Although this parameter may affect 
the nonlinear saturation of the tearing mode, the linear 
growth rate was shown to be independent of it^® . 

The choice of equal temperatures for ions and electrons 
implies that their drift speeds have the same magnitude 
but opposite directions: Uoi = —Ujje = Uuz. This also 
implies^® that the electron contribution dominates the 
growth of the tearing mode instability. The frequency 
ratio ujpf,jD,ce and the halfwidth L are related with the 
electron thermal speed Vth,e = and the initial 

drift speed, respectively, by means of 


iOp^ 1 C j ^ 

—^ =-and — =-, 

^ce 2 Vth,e L/ 2 Vth,i 


( 5 ) 


with Pi = Vth,i/^ci being the thermal ion gyroradius in 
the asymptotic magnetic field B^oy In the results to be 
shown, lengths are normalized to the ion inertial length 
di = c/ujpi and times to 

For the simulations, we use the PIC code ACRONYM 
(see Ref. 30 for a description), which allows to select dif¬ 
ferent shape functions, including all those described in 
the Appendix A. We utilize the code in its 2D3V mode, 
i.e.: no variations are calculated in the translational in¬ 
variant z direction parallel to the current. The boundary 
conditions for particles and electromagnetic fields are re¬ 
flecting in the ice direction and periodic in the ±y direc¬ 
tion. 

Other numerical parameters used in the simulations 
are the following. The grid size Ax is one Debye length 
Xoe = \/ eofcs7e/(?T-oe^) (in terms of our length unit, 
(c/wpi)/Ax = 76.8). The time step At = 0.087a;“g^ is 
chosen small enough to solve the electron motion and 
to fulfill the Courant condition AxjicAt) = 0.5 < 1. 
Following Ref. 28, the simulation box is x Ly = [1048 
Xoe] X [2048 A£)e[=[13.3 c/iOpi] x [26.6 c/ujpi], i.e.: the 
number of grid points is 1048 x 2048. The total number of 
macroparticles per species is 7.28-10®. This corresponds 
to 40 macroparticles per cell in the center of the CS. 

The size in the y direction was chosen to allow, accord¬ 
ing to the linear theory®, the development of up to seven 
unstable tearing modes with wavelength A = Ly/M, with 
M an integer satisfying 2TTML/Ly = kyL < 1 (for our 


parameters, 2tt LjLy = 0.136). Thus, we can investigate 
the interaction and exchange of energy between magnetic 
islands of different size, i.e.: multimode tearing. Note 
that for a simulation box size that would allow only one 
tearing island, this would reach stabilization and satu¬ 
ration stage at very small amplitudes, unimportant for 
space plasmas®^. 


III. ORIGIN AND CONSEQUENCES OF THE 
NUMERICALLY GENERATED ANISOTROPIC HEATING 

A. Numerical heating and shape function 

In this section, we show the results of five simulation 
runs with different combinations of shape functions (CIC 
or TSC schemes) and macroparticle number (ppc). All 
the other parameters are the same. They are denoted as 
CIC-40ppc, CIC-160ppc, CIC-360ppc, TSC-40ppc and 
TSC-160ppc. The run CIC-40ppc uses 40 macroparticles 
per cell in the center of the CS, as previously described, 
and the other CIC-160ppc and CIC-360ppc use 4 and 9 
times more particles, respectively. 

The ACRONYM code is based on the momentum¬ 
preserving scheme for PIC codes. Hence, due to the 
discretization error and the distribution of macroparti¬ 
cles according to the shape function, the total energy 
Et is almost but not completely conserved^®. As men¬ 
tioned before, the small deviation is due to the unphysi¬ 
cal non-conservative forces experienced by the macropar¬ 
ticles (see 7.6 of Ref. 20). For this reason, our first check 
is the energy conservation. In PIC codes, this quantity 
is computed as follows. 



The terms in the right hand side are the magnetic, elec¬ 
tric, kinetic and thermal energy, respectively. The sum 
is over all the species s. 14 = (l/?^s) / vfgdv^ is the 
bulk velocity, Ug = f fg dtp the number density, fg the 
distribution function and mg the mass of each species. 

The results of our simulation runs for the evolution 
of the total energy are depicted in Fig. 2(a). It can be 
seen that beyond t > 40 the conservation of en¬ 
ergy is improved by two orders of magnitude using the 
TSC over the CIC shape functions for the cases with 
40 macroparticles per cell: the total energy for TSC- 
40ppc increases above its initial value by about 0.25%, 
while for CIC-40ppc by about 25%. A similar effect is 
seen for the case with 160 macroparticles per cell: the 
total energy for TSC-160ppc increases above its initial 
value by about 0.05%, while for CIC-160ppc by about 
8%. Fig. 2(a) also demonstrates that an increase of the 
number of macroparticles does not have such strong ef¬ 
fect, while slows down the simulation to a large extent. 
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In fact, since the noise level scales as the square root of 
the number of macroparticles a reduction of numerical 
noise to 50% requires four times the number of macropar¬ 
ticles and thus it is four times slower (see e.g, CIC-160ppc 
compared with CIC-40ppc). However, the choice of TSC 
over CIC costs only a fraction of that effort and it slows 
the run by only ~ 32%. 





tiici 

Figure 2. Time history of some quantities for the runs with 
different combinations of shape functions (CIC and TSC) and 
number of macroparticles per cell (40-360 ppc). It also in¬ 
cludes a run with the XOOPIC code (using CIC shape func¬ 
tion and 40ppc) for validation purposes (denoted by A). The 
plots show: a) total energy (normalized to its initial value), 
b) average electron temperature Te,x, c) average electron tem¬ 
perature Ts,y and d) average electron temperature Ts,z- The 
temperatures are normalized to their equilibrium (initial) val¬ 
ues. Note that we break the plots in the vertical direction in 
two parts, with larger scales in the upper parts. This is in 
order to visualize more easily the small variations of the quan¬ 
tities obtained using the TSC shape function, compared with 
the large variations in those obtained using the CIC shape 
function. 

The apparent bad energy conservation of the runs with 
the CIC scheme deserves an additional validation test. 
For this sake, we also compare with the results obtained 
by using the XOOPIC code (Ref. 32), that only provides 
linear weighting (CIC). The parameters of this run are 
identical to the ACRONYM run CIC-40ppc. One can 
see in Fig. 2(a), the total energy evolution obtained by 
ACRONYM and XOOPIC. They are similar in the be¬ 
ginning, and start to diverge after trid ^ 16. From that 
time onwards, XOOPIC has a better performance (the 
energy is conserved in < 8% for the last time depicted 
tOci ~ 33) than ACRONYM if the CIC scheme is used 
(energy conservation < 16% for tfld ~ 33). We checked 
that all the physical and numerical processes to be dis¬ 
cussed in the next sections develop almost identically for 
either of the two codes, at least until the aforementioned 
time. This allows us to conclude that the ACRONYM 
runs with the CIC shape function predict a similar be¬ 
haviour of the modeled system compared to other stan¬ 
dard PIC codes. Note that we have not obtained results 


with the XOOPIC code for times greater than tfld 33, 
since we were using the free serial version of the code, 
which is much slower than ACRONYM. 

Although not shown here, we also checked that the 
use of even higher order shape functions, such as PQS, 
does not improve significantly the conservation of energy 
in comparison with TSC (while slowing down the runs 
by ^ 71% compared to the CIC scheme). This is in 
agreement with another recent study (Ref. 21). Those 
authors found that the influence of nonphysical processes 
due to the choice of the shape functions are much less 
notorious comparing cubic and quadratic interpolations 
rather than quadratic and linear ones. 

We also found that the insufficient conservation of en¬ 
ergy in the simulations using CIC shape functions is 
mainly due to artificial electron heating. As expected 
due to their larger inertia, the numerical ion heating is 
negligible. This can be seen in Fig. 2(b), depicting the 
evolution of components of the electron temperature Te.x 
averaged over the whole simulation box. The evolution 
of the electron temperature in the other two directions is 
displayed in Figs. 2(c) and 2(d). All those curves follow 
the same trend as the total energy: the TSC-runs keep 
the temperatures more constant than the CIC-runs, while 
an increase in the macroparticle number does not help to 
the same extent. Hence, we notice that a higher order 
shape function allows to resolve better the details of the 
electron motion than a simple brute force increase of the 
macroparticle number. 

It is important to note that this numerical electron 
heating is anisotropic: the increase in the out-of-plane 
temperature Tg.z is less notorious than Te^x or Tg^y (cf. 
Figs. 2(b)-2(d)). We will explore the origin of this be¬ 
haviour in the next subsections. 

Therefore, we showed that the choice of a higher order 
shape functions leads to a better conservation of energy 
and reduction of numerical heating in a more efficient 
and computationally less expensive way than the usu¬ 
ally used increase of number of macroparticles per cell. 
This reduction of numerical heating is in agreement with 
the findings of Cormier-MicheR^ for simulations of laser 
Wakefield accelerators. Note that the role of shape func¬ 
tions was already discussed in early electrostatic PIC sim¬ 
ulations comparing NGP and CIC schemes (see Ref. 33). 


B. Suppression of tearing instability by bifurcation due to 
anisotropic heating: "initially isotropic runs". 

Here we discuss the physical consequences on the CS 
evolution of the anisotropic numerical heating seen in the 
runs with CIC shape function and/or an insufficiently 
small number of macroparticles per cell. Previous stud¬ 
ies have shown that a temperature anisotropy can gen¬ 
erate bifurcated Temperature anisotropies, 

generated e.g. in simulations with CIC shape functions, 
can produce bifurcation in the out-of-plane component of 
the total current density A comparison for a spe- 
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cific time between the runs CIC-40ppc and TSC-40ppc is 
shown in Fig. 3. Jz shows a double peak structure in the 
first run but not in the second one. Instead of that, TSC- 
40ppc shows several magnetic islands. This represents a 
growth of the tearing mode. 
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Figure 3. Contours of the out-of-plane component Jz of the 
total current density, at the time t = 33 , for two runs 

a) CIC-40ppc and b) TSC-40ppc. Black lines are the mag¬ 
netic field lines. The run CIC-40ppc shows a bifurcated struc¬ 
ture while TSC-40ppc shows only a single peak structure with 
large tearing mode islands. 


The different behaviour of the evolution of CS with 
CIC and TSC shape functions can also be shown by 
means of the evolution of reconnected flux T = J i? • dS. 
In a 2D setup, this quantity can be obtained as the dif¬ 
ference of the vector potentials Az{x = 0) between the 
O and X points of each magnetic island. Since the num¬ 
ber of magnetic islands varies with time, and at the end 
there is only one left, we chose to compute this quan¬ 
tity by taking the difference between the maximum and 
minimum value of Az along the line a; = 0. This is a rep¬ 
resentative quantity when the islands can be easily dis¬ 
tinguished, although not accurate at the beginning when 
there are many small magnetic islands arising from nu¬ 
merical noise. The evolution of the reconnected flux 41 
for those cases is displayed in Fig. 4. 
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The fast growth of the reconnected flux in Fig. 4(b) 
after t ^ 60 is the signature of the explosive phase 
of reconnection, after which there is only one remaining 
magnetic island in the whole simulation box. Afterwards, 
a saturated stage is reached where the entire structure 
of the CS is disrupted due to the outflows in opposite 
directions. Note that the latter is a numerical artifact due 
to the limited size of the simulation box and the periodic 
boundary conditions in the y direction. The values of 
the maximum reconnected flux d' are of the same order 
of magnitude as reported in previous studies (see, e.g.. 
Ref. 26). Note that the reconnected flux in simulations 
with CIC shape function does not increase significantly. 
Instead, it always remains around the same value even 
for very long times, no magnetic islands are formed at 
the center of the CS, and therefore there are no X nor O 
points. This issue will discussed in more detail in Sec. V. 

In order to see more clearly the evolution of this bifur¬ 
cated structure, we plot the integrated profile of Jz across 
the inhomogeneous x direction for the CIC-40ppc run in 
Fig. 5. Bifurcation starts at around t ~ 15 As it 
was shown in Ref. 28 and confirmed by our simulations, 
it is mainly driven by a reduction of the electron current 
dentiy Je^z- This in turn, it is due to a reduction of the 
electron bulk velocity 14,while the electron density dis¬ 
tribution TT-oe does not vary too much (Je,z = —enoeVe,z)- 
The ion current density Ji z is practically unchanged. 



X/(C/C0pi) 


Figure 5. Evolution of the total current density profile 
Jz showing bifurcation (CIC-40ppc run). The profiles are 
obtained by integrating the current density along the CS: 
Jz{x) = {l/Ly) J^ll)^Jz{x,y)dy. 


Neglecting the heat flux, this is consistent with a rela¬ 
tion between the pressure anisotropy (proportional to the 
temperature anisotropy) and the derivative of the profile 
of the electron bulk velocity 14,2 (see Ref. 37) 


1 dVe^z ^ Te^x - Te,z 
By dx Tg 3 , 


( 7 ) 


Figure 4. Comparison of reconnected flux Ik for the runs a) 
CIC-40ppc and b) TSC-40ppc. The left panel a) (with CIC 
shape function) correspond to the case with anisotropic heat¬ 
ing and CS bifurcation, while the right panel b) (with TSC 
shape function) correspond to a tearing instability without 
significant anisotropic heating nor bifurcation. 


This relation also agrees with the behaviour seen in the 
TSC-40ppc run: since there is no significant electron 
heating, no significant reduction in Ve^z is observed at 
all and, therefore, no bifurcation. And, reconnection, 
the last stage of the tearing island merging, occurs for 
TSC-40ppc but not for CIC-40ppc. 
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Let us now discuss the specific effects of the tempera¬ 
ture anisotropy on the instabilities of the CS. The elec¬ 
tron temperature anisotropy can be quantified as Ag = 
T’e.ll/T’e.J-, with Tejl = Te,y and Te,± = (Te.x + Te^z)l2 
the temperatures in the directions parallel and perpen¬ 
dicular to the asymptotic Harris magnetic field (in y di¬ 
rection). Note that the temperature that enters in the 
Eq. (3) (defining the Harris CS equilibrium) is the per¬ 
pendicular one, Te^±. 

Fig. 6 depicts the evolution of the electron tempera¬ 
ture anisotropy A^, for the five “initially isotropic” runs. 
Note that the anisotropic heating (enhanced by using a 
CIC shape function) develops in the “right" physical di¬ 
rection: the parallel temperature increases above the 
perpendicular T± during the evolution of a Harris CS. 
This stabilizes the tearing mode^^. The reason for the 
preferential parallel heating can be easily understood: 
since the tearing mode develops an electric field parallel 
to the magnetic field, the electrons can be heated along 
that direction^®’^®. The physical anisotropic heating for 
the TSC-40ppc run, however, is much smaller than the 
numerical heating for CIC-40ppc. 



Figure 6. Time history of electron temperatures anisotropies 
Ae = Teji/Te^i for the five (initially isotropic) runs with dif¬ 
ferent shape functions and number of macroparticles per cell. 
The theoretical thresholds for tearing mode stabilization are 
depicted by -|- Eq.(9) and Eq.(8). The latter one is calcu¬ 
lated with the instantaneous value of Ti.u/Te.x for the run 
CIC-40ppc. 


All the runs which developed sufficiently strong tem¬ 
perature anisotropy led to bifurcation and stabilization 
of tearing mode. The reason is the threshold for the 
tearing mode stabilization in dependence on the electron 
temperature anisotropy^®. This condition can be written 


Ae 


T, 


e.ll 


T, 


e,± 


> 




Te.± 


-1 


( 8 ) 


For Te,± 3> Ti^±, it is possible to recover approximately 
the known stability criterion originally derived by LavaH® 
and Forslund^® 


Ae 


Te 


> 


T, 


e,± 



( 9 ) 


For our equilibrium parameters with x = Te^x, the 
right hand side of Eq. (8) is 1.159. But, due to the strong 
electron heating, the factor (1 -|- Ti^±/Te^±) deviates con¬ 
siderably from its equilibrium value (equal to 2) in the 
course of the CIC shape function simulations. For that 
reason, we plot Eq. (8) in Fig. 6 on dependence of the 
instantaneous value of (1 + Ti^±/Te^±) for the run with 
strongest heating CIC-40ppc (this is the minimum value 
for the threshold), instead of its equilibrium value. As 
it can be seen in Fig. 6, there is a good match with the 
stabilization threshold for all cases that developed bifur¬ 
cation and suppressed the tearing mode. 

The theoretical reason behind this stabilization (see 
Refs. 38 and 42) is the apparition of another instabil¬ 
ity driven by electron thermal anisotropies: the Weibel 
unstable mode®^’^®. According to the classical theory 
for a homogoneous plasma^^, the Weibel mode is a very 
low frequency electromagnetic wave instability (real fre¬ 
quency w ^ 0). The propagation direction and maximum 
growth rates of the Weibel instability are perpendicular 
to the warmer temperature (associated to Tg^x in our 
case), while it is damped and does not propagate in the 
colder direction (associated to Tg ||. In our 2D case, the 
direction of Weibel damped mode turns out to be the 
same as the wave vector ky of the tearing mode (y or 
parallel ||). Both Weibel and tearing modes are driven 
by Landau resonance. Hence, they are coupled and as a 
result, the tearing mode instability becomes damped as 
well. 

Note that our choice of parameters is especially sensi¬ 
tive to temperature anisotropies, since they are enhanced 
most for small gyroradii and sufficiently realistic (large) 
mass ratios®®. It is also important to remark that temper¬ 
ature anisotropy effects are relevant mostly for electrons 
at ion time scales, comparable with Moreover, the 

Weibel instability driven by ion temperature anisotropy 
is much weaker than that driven by electron anisotropy. 
Therefore, in order to determine the threshold of tearing 
mode stabilization, it is not necessary to take into ac¬ 
count ion temperature anisotropies unless they are very 
large and the electron temperature anisotropy is negli¬ 
gible smalR® . Since this is not the case in any of our 
simulations, the ion temperature anisotropy can be safely 
neglected for all the cases of simulations presented in this 
paper. 


C. Numerical CS bifurcation and entropy 

In the course of rising anisotropic heating and CS bifur¬ 
cation, a previous study (Ref. 28) analyzed the entropy 
growth. By using the XOOPIC code®^, that provides 
only the CIC shape function, it was found a growing 
entropy since the very beginning and its associated bi¬ 
furcation later on. This was associated with the phys¬ 
ical process of electron chaotic scattering which takes 
place if the heated electrons cross the center of the 
(^g46.47 .^g^g beginning of this section, 
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our run CIC-40ppc with the ACRONYM code repro¬ 
duces this bifurcation. Regarding the entropy, in Ref. 28 
the relative entropy with respect to the initial velocity 
distribution function, sometimes called Kullback-Leibler 
divergence^®, was used. Let us use instead the conven¬ 
tional definition of information entropy"^®, 

S{t) = - j f{v, t) \nf{v, t) (fv, (10) 

which differs from the thermodynamic (or Gibbs) entropy 
only by the Boltzmann’s constant kg. This quantity is 
approximated by means of a histogram estimator^® 

71a; )71y jTla: 

S' « - ^ pij^k Inpzj.fe -I- 31n(Au), (11) 

where Pij^k = fi,j,k^v^ are the bin-related probabilities 
of the macroparticles to occupy a cell (labelled i,j, k) in 
the velocity space Vx,Vy,Ve^z with a bin size Au®. In (11), 
nx,ny,nz are the number of cells in each direction of the 
velocity space and fij^k is the discrete approximation of 
the continuous distribution function f{v,t) in each cell 
labelled i-j-k. The entropy (11) is defined up to a con¬ 
stant offset depending on the relative units of Av (in our 
case, in units of the speed of light c). Since we are using 
natural logarithms, the units of S are nats (1.44 bits). 
We checked that S is more or less invariant through a 
wide range of choices of Av. Thus, the time evolution of 
the entropy for the five isotropic runs is shown in Fig. 7. 
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Figure 7. Time histories of electron (a) and ion (b) entropies 
S{t) calculated according to Eq. (11) for the five isotropic 
runs. This entropy can quantify the numerical collisions. The 
electron entropy curves follow the same trend as the curves for 
the total energy and electron temperatures shown in Fig. 2. 
And for the same reasons explained in that Fig. 2, we break 
the y axis in two halves with different scales. 


Note that the ion entropy varies very little for different 
shape functions and number of macroparticles, since the 
ion time scales are much larger than those of the elec¬ 
trons. Hence, the ion contribution to the total entropy 
is negligible. On the other hand, the electron entropy 
for the CIC-40ppc run shows a growing behaviour, re¬ 
producing the result of Ref. 28. But choosing a higher 
order scheme for the shape function (see curve for the 
TSC-40ppc run) keeps the entropy more constant in time. 
To a lesser extent, the same effect is achieved by using 
more macroparticles per cell. This is in agreement with 


our previous results concerning the suppression of quasi- 
collisional heating by using higher order shape functions 
(see Fig. 2). 

A collisionless Vlasov plasma by definition describes 
reversible processes and thus keeps the entropy constant 
for a closed system (see Ref. 51 for a discussion about 
non-equilibrium processes and the validity of the Vlasov 
approximation in a real plasma). But an increase in the 
entropy is expected in PIC codes, because the discretiza¬ 
tion of space due to the introduction of a grid and fi¬ 
nite time stepping (the coarse-graining of phase-space®^) 
can create entropy. This is true even for a system in 
thermodynamic equilibrium which should have maximum 
entropy for some given constraints (see section 12.6 of 
Ref. 19). The rate of entropy increase can be calculated 
from the Boltzmann equation with an effective collision 
operator given by Eq. (Bl). The latter is proportional to 
the increase of kinetic energy according to 


d 1 
dt 2 


(F-F)2 = (12) 

^— f dvdv'{vi-Vi)Qij{vj-Vj)f{v)f{v'). 

th,e ^ 


2v. 


This is valid for an instantaneous Maxwellian distribution 
function. Here, Qij is a tensor proportional to the right 
hand side of the numerical collision operator, therefore 
having the same dependence as Eq. (Bl). And to the 
order and approximations used in the kinetic equation, 
the rate of change of kinetic energy is equal to the rate 
of change of total energy. This explains the correlation 
between the curves for the total energy. Fig. 10(a), and 
the electron entropy, Fig. 7(a) (in particular the growing 
curves for the CIC runs). This correlation has already 
been noticed in early uniform and electrostatic simulation 
studies (see, e.g.. Ref. 53). 

The results showed in this section can be interpreted 
in the following way. It is known that the entropy in¬ 
creases due to the presence of dissipation. This process 
can be characterized by the diffusion coefficient Dij that 
appears in the effective collision operator of the Boltz¬ 
mann equation (B2). It was mentioned that this coef¬ 
ficient depends on the shape function among other pa¬ 
rameters specific to PIC codes. Therefore, the increase 
of entropy in momentum-conserving PIC codes is an in¬ 
dicator of the strength of numerical collisions. Finally, 
those numerical collisions or scattering generate the non¬ 
physical heating observed in our simulations (the second 
numerical heating mechanism mentioned in Sec. I). Note 
that previous investigations®''^ have shown that the differ¬ 
ence between the ideal Vlasov and Boltzmann equation, 
given by the numerical collision operator, is especially 
important for 2D3V plasma models in comparison with 
the 3D cases. 

The numerical collisions can even provide an expla¬ 
nation for the location in which bifurcation takes place 
(at the center of the CS). In fact, in a previous work®® 
it was shown that numerical collisions in electromagnetic 
PIC code simulations are reduced in regions with stronger 



















magnetic field. Since the magnetic field strength is min¬ 
imum, initially even zero at the center of a Harris-type 
CS, we expect that the numerical quasi-collisional effects 
are dominant there. In this way, the enhanced numer¬ 
ical collisionality produces a stronger numerical heating 
around the center a; = 0 (not shown here), which in turns 
generates a bifurcated CS. This same previous work®® no¬ 
ticed that the numerical collisions are highly anisotropic 
in a 2D configuration with an externally applied magnetic 
field. The electron-ion collisions, measured through the 
temperature relaxation time, drag and diffusion coeffi¬ 
cients, depend on the relative direction of the magnetic 
field relative to the electric field fluctuations. They stated 
that is due to the neglect of the spatial variations in the z 
direction: a 2D3V PIC code constrains the motion of the 
macroparticles to the x-y plane, but solves for all three 
components of the velocities Vx-Vy-Vz- However, the anal¬ 
ysis in Ref. 55 was carried out for a plasma embedded in 
a strong magnetic field. Therefore, it may not apply to 
our case, especially near the center of the CS. For this 
reason we tested that behaviour by means of two small 
simulations with CIC shape function to enhance the nu¬ 
merical collisions. The first one is a 2D run with the same 
parameters as the run CIC-40ppc, a smaller simulation 
box [128Aa; x 128Aa;] and without the initialization of 
a Harris sheet (only a constant magnetic field in the y 
direction with magnitude B^oy)- And the second one is 
a full 3D simulation with the same previous parameters 
but in a box [128 Ai x 128 Ai x 128 Ax]. The 3D run 
developed a much smaller numerical heating and electron 
temperature anisotropy than the 2D one. This indicates 
that 2D simulations overestimate the effects of numerical 
collisions and thus, the origin of the anisotropy is mainly 
due to the neglect of 3D effects. 


D. Numerical heating and background plasma effects 


So far, we did not include a background population 
in the previously described runs. This might have im¬ 
portant implications in the numerical scattering, since at 
the edge of the CS there are large electromagnetic fluctu¬ 
ations due to the small number of macroparticles in that 
region. In order to clarify if the neglect of a background 
plasma is affecting our conclusions, we investigated its 
influence in our simulations. 

For our standard runs with 40ppc and CIC scheme 
the results are depicted in Fig. 8. These runs are for 
the cases with no background (CIC-40ppc), 10% (4 ppc, 
CIC-40ppc-back01) and 20% (8 ppc, CIC-40ppc-back02) 
of the peak density of the Harris plasma population. 
All of them show a similar evolution for the temper¬ 
ature anisotropy. Correspondingly, we found, by look¬ 
ing the time history of Jz profiles, that the bifurcation 
of the CS evolves in a practically indistinguishable way 
compared to the case without background (see Fig. 5). 
Note, however, that the total energy is less conserved for 
higher plasma background densities. This could be inter¬ 


preted as an indication that the influence of the scatter¬ 
ing is enhanced mainly due to the electron temperature 
anisotropy, and does not depend too much on other in¬ 
volved quantities that may be affecting the evolution of 
the total energy. 
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Figure 8. Time history of a) total energy (normalized to its 
initial value) and b) electron temperature anisotropy Ae = 
Teji/Te.x for the cases of a CS without background (CIC- 
40ppc), background population at 10% of the peak density 
(CIC-40ppc-back01) and at 20% of the peak density (CIC- 
40ppc-back02). In all those runs the CIC scheme has been 
used with 40ppc. 


As an additional test, we plotted the evolution of the 
total energy and electron temperature anisotropies for 
simulations with the same parameters used in our stan¬ 
dard runs shown in Fig. 2, but with the addition of a 
background population of 10% of the peak density. The 
results are shown in Fig. 9. By comparing with the cor¬ 
responding cases without background plasma (Figs. 2(a) 
and 6), one can see an almost identical evolution for the 
curves of the CIC scheme. Only for the TSC scheme 
some differences are obtained in the evolution of the elec¬ 
tron temperature anisotropy. This is because in those 
cases, the explosive phase of magnetic reconnection (more 
accurately modeled with the second order interpolation 
scheme TSC) takes place a little bit earlier in the back¬ 
ground cases than without it. This is in agreement with 
previous Vlasov (Ref. 56) and PIC (Ref. 42) simulation 
results of magnetic reconnection. Those authors reported 
that the onset of magnetic reconnection and the satura¬ 
tion of tearing instability vary little depending on a back¬ 
ground, while the reconnection rates are substantially re¬ 
duced with an increasing background plasma density. We 
also observed that the explosive phase of reconnection is 
shorter without background than with it. This is due to 
the larger Alfven velocity in the low density regions away 
from the center of the CS, where it decreases exponen¬ 
tially. 
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Figure 9. Time history of a) total energy and b) electron 
temperature anisotropy Ae = Tg, 11 /^ 6 . _l for our standard runs 
described in Sec. Ill A, but with the addition of a background 
population of 10% of the peak density. Compare with the 
Figs. 2(a) and 6. 

Therefore, our conclusion is that the addition of a back¬ 
ground population does not seem affect too much the bi¬ 
furcation/stabilization of the tearing mode in the course 
of the evolution of a run with a CIC scheme. This imply 
that nor the developed increasing temperature anisotropy 
nor the numerical scattering are affected too much by this 
population. Hence, the large fluctuations at the edge of 
a current sheet (due to the limited number of macropar¬ 
ticles in the simulations without background plasma) do 
not affect too much the physical processes that we mod¬ 
eled in this paper. This is because bifurcation takes place 
near the center of the current sheet, away from the re¬ 
gions with large fluctuations due to the too small number 
of macroparticles per cell. For other physical processes, 
the inclusion of a background may become critical for 
a correct description (e.g.: any estimation that requires 
the calculation of momenta of distribution function at 
the edge of the CS). 


IV. STABILIZATION OF TEARING MODE FOR 
INITIALLY IMPOSED TEMPERATURE ANISOTROPIES: 
"INITIALLY ANISOTROPIC RUNS" 

A. Initial temperature anisotropy relaxation 

Let us investigate the cause for the relation between 
anisotropic heating and CS bifurcation. For this sake, 
we performed simulations with an initial temperature 
anisotropy A^,, obtained by choosing T^ y = Tg || below 
and above its equilibrium value Te^±. All the other pa¬ 
rameters are kept identical to those of the TSC-40ppc 
run, in order to minimize the numerical heating effects, 
in particular = Te,x = Te^z- Note that due to 

our choice of Ti = Tg, we also have an ion temperature 
anisotropy Ai = These additional five runs, 

with temperature anisotropies A^ = Ai = 0.64, 1.21, 
1.44, 1.69 and 1.96 (corresponding to the choice of elec¬ 
tron thermal speeds of Vth,ey/vth,ei. = 0.8, 1.1, 1.2, 1.3 
and 1.4, respectively), are denoted as TSC-40ppc-A0.64, 
TSC-40ppc-A1.21, TSC-40ppc-A1.44, TSC-40ppc-A1.69 
and TSC-40ppc-A1.96, respectively. We preferred to 
sample more values with Ag > 1 than Ag < 1, be¬ 


cause those correspond to the temperature anisotropies 
obtained using CIC shape functions (cf. Fig. 6). 

Some results of the ’’anisotropic runs" are illustrated 
in Fig. 10. The parallel temperature Tg ^, shown in 
Fig. 10(b), tends to decrease towards its equilibrium 
isotropic value for all the runs with Ag > 1. The fi¬ 
nal value is proportional to the initially imposed Ag. On 
the other hand, we can see in Fig. 10(a) that Te^x (same 
for Tg^^, not shown here) increases departing from its 
initially imposed equilibrium value, approaching to an 
asymptotic value also proportional to the initial value of 
Ag. In summary, in all those cases the system relaxes 
towards a more isotropic state during a short transient 
time, as can bee seen in Fig. 10(c). 
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Figure 10. Time history of quantities characterizing the 
five runs with initially imposed anisotropy Ag (’’anisotropic 
runs”), compared to the “isotropic run” TSC-40ppc. a) av¬ 
erage electron temperatures b) average electron tem¬ 

peratures = Tg^ii and c) average electron temperature 
anisotropy Tg^i/Tg^^ for different initially imposed tempera¬ 
ture anisotropies Ag G [0.6, 2.0]. The theoretical thresholds 
for stabilization of tearing mode are shown with lines -|- Eq.(9) 
and * Eq.(8). The latter one is calculated with the instanta¬ 
neous value of Ti,x/rg,x for the run TSC-40ppc-A1.96. 

At a first sight, the results regarding the fast relax¬ 
ation of the initially imposed anisotropy shown in Fig. 10 
may seem to be a numerical artifact. This is because, as 
we already mentioned, the only temperature that enters 
into the Harris sheet equilibrium is the perpendicular one 
Tg^x- Therefore, an initial state with a different paral¬ 
lel temperature Tg y is still a Harris equilibrium, and the 
initially imposed temperature anisotropy should be con¬ 
served in an ideal collisionless Vlasov system without ad¬ 
ditional instabilities (although this equilibrium turns out 
to be unstable to an anisotropy driven instability under 
the conditions to be described). Hence, in order to prove 
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the validity of the runs shown in Fig. 10, we carried out 
two simulation runs for the case with the highest consid¬ 
ered anisotropy (TSC-40ppc-A1.96), while changing the 
number of macroparticles per cell: TSC-360ppc-A1.96 
(360 macroparticles per cell) and TSC-lOOOppc-Al.96 
(1000 macroparticles per cell). Those runs have, respec¬ 
tively, three and five times less numerical noise compared 
to TSC-360ppc-A1.96, since (as we already mentioned) 
the noise level is reduced according to y/N, with N the 
number of macroparticles per cell. 


H 
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Figure 11. Time history of the average electron temperature 
anisotropy As = Ts^\\/Ts^± for cases with an initially imposed 
temperature anisotropy As = 1.96, TSC shape function and 
different number of macroparticles per cell (40, 360 and 1000). 
For comparison, the result for an isotropic initial distribution 
{As = 1.0) is also plotted with a solid line. 


One can see in Fig. 11 that in case of more macropar¬ 
ticles per cell, the electron temperature anisotropy Ag 
does not fall as quickly as for the case TSC-40ppc-A1.96, 
especially for later times. There is a trend (convergence) 
towards the reduction of the relaxation of temperature 
anisotropy for an increasing number of macroparticles per 
cell (see, e.g., the curve for the run TSC-lOOOppc-Al.96 
in Fig. 11). This behaviour can be attributed to the nu¬ 
merical collisions, since they enhance the isotropization of 
the distribution function (even by using the TSC scheme) 
and are weaker for a larger number of macroparticles per 
cell. This process is developed via pitch angle scattering, 
redistributing the particle velocities in the momentum 
space: a momentum transfer takes place from the par¬ 
allel velocity direction to the perpendicular one, with a 
negligible change of the particle energy (total speed). 


B. Physical origin of the initial temperature anisotropy 
relaxation: Weibel instability 

In spite of the previous numerical argument, the initial 
drop in Ag = Tgji/Te^j, shown in the curves of Fig. 10(c) 
or Fig. 11 (especially for high values of Ae) seems to be 
independent on the number of macroparticles per cell, in¬ 
dication of a physical origin. In order to understand this 
process, it is convenient to look more closely the evolu¬ 
tion of the electron temperature anisotropy in this initial 


stage, focusing in the case TSC-360ppc-A1.96. This is 
shown in Fig. 12 
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Figure 12. Time history of the out of plane magnetic field 
and temperature anisotropy averaged along y, at different dis¬ 
tances X from the center of the CS, for the run TSC-360ppc- 
A1.96. a) and b) Magnetic field \Bz\. c) and d) Electron 
temperature anisotropy Ts^y/Ts^x- The left plots a) and c) 
are for the early stages of the system while the right plots b) 
and d) are for the long term evolution. 


We can distinguish three stages in the evolution of the 
system: (1) the initial between tVld = 0 — 0.5 in which 
the anisotropy is constant, (2) between tVld = 0.5 — 1.5 
in which the anisotropy is quickly relaxed and (3) t^ld > 
1.5, where the system relaxes its anisotropy slowly until 
very late times. 

In the initial stage (1), due to the high value of the 
initial imposed anisotropy Te^y > Te^x,Te^z, the unmag¬ 
netized center of the CS becomes unstable to the Weibel 
instability. The timescales in which Weibel operates are 
much shorter than these of the tearing mode, and there¬ 
fore the aforementioned coupling between these instabil¬ 
ities can be neglected in this early stage (as well as any 
consequence of this coupling, such as the stabilization 
of tearing mode). We already mentioned that this in¬ 
stability propagates in the colder direction associated to 
Te^j_ = Ts^x (since no wave vectors in z direction are al¬ 
lowed in our 2D configuration). Its growth rate can be 
approximated as 



This expression is valid in the low frequency regime 
{uj/{kxVth,e) "C 1). We have also neglected the ion contri¬ 
bution, since is negligible by a factor of mi/mg. There¬ 
fore, this is a instability mainly driven by electrons. From 
the previous expression it is possible to determinate the 
threshold: 


-1 > 


kxC 


UJ. 


pe 


(14) 
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Although this electromagnetic instability was originally 
discovered and described for an unmagnetized plasma, it 
can also be triggered under the influence of a background 
magnetic held (in this context, it is sometimes called “hl- 
amentation instability”®^’^®), as long as its magnitude is 
not too high®®. Indeed, this kinetic instability only de¬ 
pends on the bulk properties of the plasma (anisotropy), 
and operates under a wide range of conditions due to 
the fact that does not depend on wave-particles reso¬ 
nances effects®®. The Weibel instability generates sta¬ 
tionary magnetic helds (which might be important, e.g., 
for the reconnection process®^) perpendicular to the di¬ 
rection of the “warmer” temperature, (in our 2D setup, 
the z direction) by the well known mechanism described 
in Ref. 62 and 63. The addition of a strong background 
magnetic held (in our case, in the y direction) tends to 
inhibit this instability because the motion of the particles 
perpendicular to it is suppressed, in such a way that it is 
not possible the momentum transfer between the paral¬ 
lel and perpendicular directions. Under the assumption 
Dee ojpe (valid in the weak magnetic helds close the 
center of the CS), it is known that the Weibel growth 
rates given by Eq. (13) are not modihed to a good ap¬ 
proximation. The background magnetic held only has 
an inhuence in the apparition of a small real frequency 
w ~ Dee (see Refs. 64-66 for further details). 

Therefore, the generation of the magnetic held (see 
Fig. 12(a)) close to the center of the CS, as well as 
the corresponding decrease in the electron temperature 
anisotropy (see Fig. 12(c)) are the main signatures of 
the development of the Weibel mode. For a quantitative 
comparison, we need a value for the wavevector kx in the 
threshold condition Eq. (14). The kx^max corresponding 
to the maximum growth rate given by Eq. (14) cannot be 
applied directly, since it may predict longer wavelengths 
than the ones allowed by the unmagnetized condition of 
the Weibel instability Dee "C ^pe- Indeed, according with 
Eq. (13), kx,max can be calculated as"^'*: 


A ^ 2L ^ Idi (to be justihed a posteriori), we can esti¬ 
mate the right hand side of Eq. (14) as (kxc/ujpe) ~ 0.21. 
This value matches well with the observed saturation 
anisotropy after tVld > 1.5. Note that this value for 
the saturation of Weibel mode is higher than the one re¬ 
quired for the stabilization of tearing mode (compare to 
Eqs. (8) and (9)). This means that a system unstable to 
Weibel mode will already be stable to tearing mode, in¬ 
hibiting reconnection. And viceversa, a PIC simulation 
that develops high enough numerically generated tem¬ 
perature anisotropy along all the center of the CS will 
stabilize tearing mode before of reaching a state unstable 
to Weibel mode (this is what happens in the previously 
discussed CIC runs). 

It is possible to validate the previous estimation for the 
saturation wavelength of the Weibel mode by looking at 
the spatial structure of B^, displayed in Fig. 13. 



X / (c/COpi) 

Figure 13. Contours of the out of plane magnetic field Bz 
for the run TSC-360ppc-A1.96. The snapshots correspond to 
three characteristic times during the evolution of the Weibel 
instability 
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with the corresponding maximum growth rate 
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We can note that the associated Xx,max = l^x,max gets 
larger for smaller anisotropies. For our parameters, the 
initial value Ae = 1-96 corresponds to Ax,max = 0.8di, 
but the value at saturation ^ 1.2 corresponds to 
Ax,max = l.Sdi. This length would reach strongly magne¬ 
tized regions, suppressing the Weibel mode. In Ref. 67, it 
was proposed that Weibel instability will saturate when 
the wavelength along x is of the order of magnitude 
of the halfwidth L of the CS (measured from the cen¬ 
ter), since it is the typical lengthscale in which the mag¬ 
netic held varies. Assuming that reasonable wavelength 


The “chessboard“ (hlamented) structure around the 
center of the CS has already been observed in previ¬ 
ous studies of magnetic reconnection, in both electron- 
positron plasmas (Refs. 68-70) as well as in electron- 
proton plasmas (Refs. 67 and 71). In those works, it 
has been identihed as a characteristic signature of the 
Weibel mode generated by the temperature anisotropy 
inside of the magnetic tearing islands. In this con¬ 
text, the anisotropy is caused by the particles bounc¬ 
ing inside of those structures, i.e.:, the well known Fermi 
mechanism^^. The Weibel instability has also been pro¬ 
posed to be generated inside of the electron inertial length 
around the X point, as a result of the electron inhow"'’®’®^. 

At the beginning, when the anisotropy is high, the size 
of the Bz structures across x is relatively small, as can 
be seen in Fig. 13(a). They reach larger sizes later, when 
the anisotropy is reduced, in agreement with Eq. (15) 
(see Fig. 13(b)). Their saturation size at tfld ~ 1.5 is 
of the order of A = Id^, as expected from our previous 
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estimation. They cannot grow further because of the 
constrain in the magnetic field strength. 

Regarding the size of the Bz structures along y, the 
usual linear theory for a homogeneous plasma cannot 
predict their length scale (there should be no wavevec- 
tor along the warmer y direction). But a mechanism 
proposed in Ref. 70 can explain and predict the finite ky 
observed in Fig. 13. A wavevector along the direction of 
higher temperature can be produced with a wavelength 
of spatial scales that matches the double of the electron 
gyroradius based on the Bz generated around the neu¬ 
tral point (see also Ref. 71). At the saturation time of 
the Weibel instability, the average component of the mag¬ 
netic held is about Bz ^ Q.\Baoy^ implying a gyroradius 
Pe,zldi ~ 0.35. This is about 1/4 of the observed size of 
the the Weibel structures at that stage (1.3(ii), consistent 
with the predictions of Ref. 70. This mechanism also ex¬ 
plains why for later times (see Fig. 13 for tQ.ci = 3.7) the 
hlamented Weibel structures in y direction are smaller 
than for tftd = 1.5, since the generated magnetic held Bz 
becomes weaker after the relaxation of the temperature 
anisotropy. This implies a smaller electron gyroradius on 
this component and thus a smaller size of the wave vector 
along the y direction. 


C. Consequence of the initial temperature anisotropy 
relaxation: Mirror instability 

After tflci > 1.5, the magnetic helds generated by the 
Weibel instability start to slowly be dissipated, as can be 
seen in Figs. 13(c), 12(a) and 12(b). This in an indica¬ 
tion of the turning off of this instability, since the Weibel 
mode generates turbulent magnetic helds®® that scatters 
the particles in a such a way that those tend to reduce 
the anisotropy in the distribution function (the source of 
their free energy). For later times, the anisotropy at the 
center of the CS a: = 0 tends to be kept at the saturation 
values, although with a slow decrease (see Fig. 12(d)). 
From this same hgure it is possible to see that the low 
density regions further away from the center, not affected 
in a signihcant way by the Weibel instability, will expe¬ 
rience a a faster isotropization in comparison with the 
center (due to the enhanced collisions). Note that this 
effect cannot be seen in the average of the electron tem¬ 
perature anisotropy (see Fig. 11). 

Although the Weibel instability seems to be active un¬ 
til tfld ^ 1.5, there is an important difference in the 
behaviour of the system after tild ~ 0.5. From that par¬ 
ticular time onwards, a secondary instability driven by 
temperature anisotropy starts to develop around the neu¬ 
tral point in the CS: the so called mirror instability 
Different from the unmagnetized Weibel instability, this 
electromagnetic instability requires a magnetic field in or¬ 
der to exist. This is precisely the magnetic field provided 
by the Weibel instability. The mirror instability is one 
of the instabilities that can be excited in a plasma when 
Tj^± > Tj ||, for both ions (classical hydromagnetic ion 


mirror mode^^) and electrons^^. Its threshold condition 
can be written as’’®^^^: 
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where the plasma betas are /3jj| = ^md I3j^j_ = 

bV(^/) • Here, the parallel (||) and perpendicular (T) 
directions are with respect to the local magnetic field. 
The well known mirror threshold condition derived from 
fluid equations^®’^® assumes cold electrons, resulting in a 
vanishing last term in (17) or (18) (this term is a correc¬ 
tion due to the field aligned electric field arising from the 
differential motion of ions and electrons^®). However, we 
checked that this last term does not play an important 
role in our this system, since initially is imposed to be 
zero due to our initialization (in addition, it saturates 
at high enough electron temperatures). It is important 
to notice that the threshold condition (17) involves the 
contribution of both ions and electrons. In our originally 
anisotropic runs, both species are anisotropic and con¬ 
tributes in more or less the same proportion to the thresh¬ 
old. The mirror modes propagates at an oblique angle 
with respect to the magnetic held (maximum growth rate 
nearly in the perpendicular direction), and have predom¬ 
inant longitudinal fluctuations^®. One key feature of the 
mirror instability is its preferential development in high 
(3 plasmas with relatively low anisotropies®®, since the 
threshold is lower (see Eq.(18)). Therefore, a weakly 
Weibel magnetized center of a CS (high fi) with an im¬ 
posed temperature anisotropy has the ideal conditions 
for the triggering of this instability®®. 

It is important to note that the (electron) whistler in¬ 
stability is other instability driven by the same condition 
Te__L > 7 / 5,11 ■ But we ruled out the existence of that in¬ 
stability in this system, since its maximum growth rate 
is in the parallel direction to the magnetic held (in this 
case, mostly in the z direction), which is not allowed in 
our 2D simulations. 

Since that around the center of the CS sheet the mag¬ 
netic topology is complex (not only a Bz is generated, 
but also a less stronger B^ component), in order to cal¬ 
culate Eq. (17) properly, it is necessary to rotate the 
tensor pressure to keep it always locally aligned with 
the magnetic held. This means that the || direction is 
imposed to be in the B/\B\ direction. In this way, we 
can see that small scale islands (with a maximum typical 
length scale of ^ c/wpe across the x direction) are formed 
around the center of the current sheet with re,-L > Te,\\ 
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(see Fig. 14), where the mirror condition Eq. (17) is ful¬ 
filled (see Fig. 15). 



- 2 . 3 - 1.2 0.0 1.2 2.3 - 2 . 3 - 1.2 0.0 1.2 2.1 

X / (c/cOpi) X / (c/cOpi) 



X / (c/cOpi) 

Figure 14. Contours of the electron temperature anisotropy 
Fe.n/re.x = Pii^e/Fx.e = /3||,e/,de.x for the run TSC-360ppc- 
A1.96. The snapshots correspond to the same three times 
shown in Fig. 13. Note that regions with Te^x > re,|| are 
shown inside of the black contour lines. 



- 0.4 - 0.2 0.0 0.2 0.4 

X/{c/(Bp|) 


X T . , 

- 1.0 1.3 3.7 6.0 

logio(D) 

Figure 15. Contours of the mirror threshold value from 
Eq. (18): logjQ(P). Only positive values are shown (corre¬ 
sponding to instability). This plot is for the run TSC-360ppc- 
A1.96 at the time tflci = 1-5, corresponding to the peak of 
activity of both Weibel and mirror instabilities. Compare 
with the corresponding plots for Bz in Fig. 13 and tempera¬ 
ture anisotropy Teji/Te.x in Fig. 15. Note that only a small 
region close to the current sheet is shown, in order to visual¬ 
ize the mirror structures that have a max;imum typical length 
scale of c/oipe across the x direction. 

The time evolution of this process can be seen in 
Fig. 16, showing the ratio of grid points unstable to the 
mirror instability in a region |a;| < L and their mean 
value, as well as the condition re,||/’?e,_L < 1 that can 
trigger this instability (number of points that satisfy this 
condition divided over the total in |a;| < L). We can 
see that the regions unstable to mirror instability start 
to appear only after tflci = 0.5. They reach their peak 
(both in size and strength) around tfld = 1.5, correlated 


with the maximum of the Weibel generated magnetic 
field Bz and the minimum in the temperature anisotropy 
Te^y/Te^x- Between those times, the source of the free 
energy of the mirror mode, the field aligned tempera¬ 
ture anisotropy Tg^x/Feji > 1, is quickly exhausted (see 
Fig. 12). The mirror instability produces strong turbu¬ 
lence in those unstable regions, deforming and kinking 
the in-plane magnetic field lines (since it acts mainly in 
the perpendicular direction to the local magnetic field). 
The mirror structures only survive for short times. Note 
that the Weibel instability by itself cannot explain the 
sudden appearance of the small scale structures seen only 
between 0.5 < tfld < 3. 




01234 01234 

t Hji t IIqi 

Figure 16. Time evolution of quantities related with the mir¬ 
ror instability for the run TSC-360ppc-A1.96. a) Number of 
grid points in |a;| < L that satisfy Tf,^]\/Te,± < 1 divided over 
the total in this region. Note that Te^n/Te.x is the electron 
temperature anisotropy aligned with the local magnetic field, 
i.e.: the || direction is in the B/\B\ direction. The choice 
of the |a;| < L is basically because encompass all the area 
that becomes unstable to mirror instability b) Mean value 
of Te.ii/Te.x for the grid points that satisfy the condition 
Te^ii/Te.x < 1- c) Number of grid points in |a;| < L that 
satisfy the mirror threshold condition Eq. (18), divided over 
the total in this region, d) Mean value of the mirror threshold 
condition (left hand side of Eq. (18)) calculated for the points 
that satisfy that condition. 


Other signatures that can support our conclusion that 
those structures are due to the mirror instability (in ad¬ 
dition to the threshold condition) can be provided by a 
comparison with the linear theory. However, due to the 
coupling with the Weibel instability acting at the same 
time, it is not easy a direct comparison, especially in our 
case of non homogeneous plasma with the mirror insta¬ 
bility only acting in a very narrow region. Nevertheless, 
the linear theory can provide us with some insight. In 
this context is meaningful to write the usual expression 
for the maximum growth rate of the mirror mode^®’^®. 
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assuming cold electrons and an ion anisotropy: 



where D is the left hand side of the mirror threshold con¬ 
dition given by Eq. (18). Because the condition of almost 
perpendicular propagation, the second term is neglected 
for the mirror mode (it gives the growth rate of the com¬ 
plementary instability: the so called firehose instability 
driven by 1 | > _l with parallel propagation). A more 
precise growth rate taking into account similar electron 
and ion contributions to the temperature and anisotropy 
of the system (more suitable for our case) was given in 
Ref. 76: 


^max — 


kLVth,i,\\Ti^\\ (l + Ti\ ) A.||)/2] 






( 21 ) 


with 



Although this expression increases with fc_L, it is only 
valid in the small Larmor radius approximation k±pi <C 
1. Following Ref. 76, the instability growth rate will de¬ 
crease in the short wavelength limit due to finite Lar¬ 
mor radius effects, and therefore would reach a maximum 
around the ion gyroradius ~ pi. In Fig. 15, the longi¬ 
tudinal size of the spatial structures (tilted with respect 
to the y direction), where the mirror condition is satis¬ 
fied, agrees well with this prediction. Indeed, a typical 
structure has a size in x ly = Z|| x l± = 0.2di x 2.0di, 
while Xj_^max = 27r/pi ^ 3di. 

On the other hand, the mirror mode predicts a small 
propagation vector along the parallel direction. Although 
most of the magnetic field has z component, there is 
a small contribution from the x component. Since our 
2D geometry does not allow wavevectors in the out-of- 
plane direction z, the parallel propagation direction has 
to be taken necessarily by the x direction (although the 
contribution for the parallel temperature comes mostly 
of Te^z)- Therefore, the mirror mode will propagate 
mostly in the perpendicular y direction in such a way 
that kj_ = ky ^ /c|| = kx- The angle 9 of propagation 
(with respect to the magnetic field) can be estimated by 
the following expression: 

tan(6») = ^ = -^. (23) 

h VA 

For typical parameters in our simulation inside of the 
mirror structures (/3jj| ~ /3e,|| ~ 80, and the field aligned 


anisotropies Ae ^ Ai ^ 0.99), this angle is around 
0 ~ 85°, which agrees very well with the tilting angle of 
a typical mirror structure {6 = atan(2.0di/0.2di) ^ 85°). 
Note that although the range of variation of the param¬ 
eters inside the magnetic mirror structures can be very 
large, especially regarding the plasma betas, the angle 
9 have a weak dependence on those parameters (mirror 
modes usually have nearly perpendicular propagation), 
and does not deviate too much from this value (different 
from the growth rate, that can be very sensitive to those 
variations). 


D. Growth rates of the initial temperature anisotropy 
relaxation due to Weibel/mirror instabilities 

In order to compare growth rates of both Weibel and 
mirror instabilities, we used the results of runs with dif¬ 
ferent levels of initially imposed anisotropy and 360ppc. 
For lower values of the initial anisotropy, the peak value 
of Bz is lower and it is reached for later times. There¬ 
fore, by performing a linear fit in the plots of the mean 
value of \Bz{x = 0)| between the times in which this 
magnetic field grows exponentially, we can measure the 
Weibel-mirror mode growth rates. The results are shown 
in Fig. 17. As we already mentioned, when the initial 
anisotropy is lower than A^ ^ 1.2, the Weibel mode is 
stabilized and it cannot generate the magnetic field nec¬ 
essary for the formation of those structures. 



n t al^y/Te,x 



n t al^y/Te,x 


Figure 17. a) Comparison among simulated growth rates 
of the Weibel instability during the linear stage on depen¬ 
dence of the initially imposed electron temperature anisotropy 
Te,ylTe,x, a fitting according to Eq. (13) (for fixed fc), a fitting 
according to Eq. (16) (for kx,m,ax), and the theoretical growth 
rates calculated from the initial parameters in Eq. (13) (for 
fixed k). See text for details, b) Saturation time tsat of the 
Weibel instability on dependence of the initially imposed elec¬ 
tron temperature anisotropy Ts^yjTs^x- 


In Fig. 17, the fitting by continuous line was calculated 
by using 'y{Ae)/Xlci = (cijAe) x (Ae - 1 - C 2 ), equal to 
the functional dependence of the analytical growth rates 
given by Eq. (13) for the Weibel instability, assuming a 
fixed k. The resulting fitting coefficients are ci = 12.3 
and C 2 = 0.20. This implies a kx corresponding to a wave¬ 
length Xx = l.Odi, in good agreement with the maximum 
size of the structures seen in the out of plane magnetic 
field B^ ■ The other fitting in dashed-point line was calcu¬ 
lated by using "f{Ae)/VLci = (ca/Ae) x (A^ - l)^/^, equal 
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to the functional dependence of the analytical growth 
rates given by Eq. (16) for the Weibel instability assum¬ 
ing a kx,max dependent on the temperature anisotropy, as 
predicted by the linear theory. The resulting fitting coef¬ 
ficient is C 3 = 10.22. It can be clearly seen that this curve 
does not reproduce well enough the behaviour of the sim¬ 
ulated growth rates. On the other hand, the predicted 
growth rates obtained by applying directly Eq. (13) for 
our parameters and the previous fixed kx (corresponding 
to the maximum size of the observed structures in B^), 
do not agree with the observed ones (see dashed line in 
Fig. 17). The last ones are about three times lower than 
those given by the linear theory Eq. (13). There are 
many reasons for this disagreement, since we are clearly 
applying this theory beyond its limits (e.g.: assumptions 
of homogeneity, weak magnetic field, combination with 
mirror instability). Therefore, we conclude that the lin¬ 
ear theory for a homogeneous plasma can only predict the 
right functional dependence on Af, for the growth of the 
Weibel instability (combined with mirror), even in our 
complex geometry of a Harris current sheet (see the fit¬ 
ting with continuous line in Fig. 17). However, the latter 
is valid only under the assumption of choosing before¬ 
hand a fixed k that depends on the typical length scale 
of variation of both magnetic field and density (double of 
the halfwidth 2L ^ 2pi ^ 1.0d^). 

On the other hand, the mirror growth rates (19) or (21) 
are very sensitive to the value of Pi , whose value can vary 
over some orders of magnitude inside of the mirror struc¬ 
tures and so is the range of variation of the predicted 7 . 
Therefore, the value of 7 of the mirror instability inside 
of the structures does not provide information useful to 
be compared, at least in an easy way, with any of the 
large scale electromagnetic fields or momenta of the dis¬ 
tribution function observed in our simulations. 

So, all the previous arguments seems to indicate that, 
in addition to the Weibel instability, a mirror instabil¬ 
ity is also acting close to the center of the current sheet 
driven by the electron temperature anisotropy. Previous 
studies have indicated the close relationship between the 
mechanism that triggers those instabilitiesEspecif- 
ically, in Ref. 81 it was noted that the linear responses of 
the distribution function for the mirror-type and Weibel- 
type perturbations, in addition to the linear growth rates, 
have the same form, replacing the quantities where the 
ion/electron gyroradius is involved (in the mirror case) 
with the electron inertial length (in the Weibel case). 

It is also is important to mention that this interplay 
between Weibel and mirror instabilities has already been 
pointed out explicitly in Ref. 82. These authors con¬ 
cluded that in a originally homogeneous unmagnetized 
but anisotropic plasma, the magnetic fields generated by 
the Weibel instability can produce mirror modes struc¬ 
tures. Those will be ordered into a chain of bubbles 
or holes, in addition to the filamented structures due 
to Weibel. The conclusion is that a homogeneous but 
anisotropic plasma will naturally evolve into a turbu¬ 
lent state due to the interaction of those instabilities (see 


Ref. 82 and references therein). 

We cannot rule out completely the presence of other 
temperature anisotropy driven instabilities in our sys¬ 
tem. Initially, the regions away from the center of the 
CS satisfy the condition to be firehose unstable. This 
is (mostly) a parallel propagation instability driven by 
the opposite condition to the mirror instability Tg y > 
But we have not found signatures of its pres¬ 
ence (e.g.: absence of resonant protons, no clear wave 
vector in those regions). The reasons seems to be that 
its predicted growth rates are about one order of mag¬ 
nitude lower than those of the Weibel mode. Besides of 
that, all the phenomena related with bifurcation happens 
close to the center of the CS, especially the reduction in 
the temperature anisotropy, where this kind of parallel 
propagating instabilities should not play a role. How¬ 
ever, it is important to remark that, under the right con¬ 
ditions, the firehose instability can also be triggered in 
a current sheet, producing kinking of the magnetic field 
lines, scattering particles and reducing the temperature 
anisotropy®'^, in a similar way to the phenomena observed 
in our case. 


E. Suppression of tearing instability by bifurcation due to 
anisotropic heating: "initially anisotropic runs". 

The main consequence of an initial electron tempera¬ 
ture anisotropy Hg > 1.2 is that, after a short transition 
period, a bifurcated CS is formed. This value corresponds 
to the one calculated in Sec. IV B for the stabilization of 
the Weibel instability, assuming the mentioned constraint 
in kx due to the inhomogeneity in the background Harris 
magnetic field. Note that this value is higher than the 
theoretical threshold for the stabilization of tearing mode 
(Eqs. ( 8 ) and (9)). Indeed, in Fig. 10(c), it is shown the 
threshold for stabilization of the tearing mode according 
to Eq. ( 8 ). Since the ratio Ti^j_/Te^± does not change 
too much away from its equilibrium values, the threshold 
curve is almost horizontal. Therefore, a system unstable 
to Weibel mode not only will be stable to the tearing 
mode (as we mentioned in Sec. IV B), but also will de¬ 
velop bifurcated structures in Jz- 

The morphology of the bifurcation in Jz is similar to 
the resulting from evolving temperature anisotropy (see 
Fig. 5). The bifurcation is stronger (i.e.: the dip in the 
current density is deeper) for higher initial temperature 
anisotropies. See Fig. 18 for a comparison of this be¬ 
haviour in the evolution of the profile Jz (x) for the runs 
with highest considered anisotropy and different number 
of macroparticles per cell: TSC-40ppc-A1.96 and TSC- 
360ppc-A1.96. Note that the initially imposed anisotropy 
is still causing a bifurcation of the CS for the cases with 
more macroparticles per cell and, correspondingly, still 
stabilizes the tearing mode instability. However, the drop 
in Jz at the center of the CS in TSC-360ppc-A1.96 is less 
deep compared with the case TSC-40ppc-A1.96. This in¬ 
dicates a trend towards the reduction of bifurcation in 
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the limit of negligible scattering. 



x/(c/copi) 
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Figure 18. Evolution of the total current density profile Jz 
showing bifurcation for runs with the highest considered ini¬ 
tially imposed temperature anisotropy Ae. = 1.96, TSC shape 
function and different number of macroparticles per cell: a) 
TSC-40ppc-A1.96 and b) TSC-360ppc-A1.96. The profiles 
are obtained by integrating the current density along the CS: 
J^{x) = (l/Ly) f^l^y^^^Jz{x,y)dy. 


In order to understand the onset of the bifurcated 
structures, it is necessary to recall the results concerning 
the relaxation of the initial anisotropy by the two insta¬ 
bilities previously discussed. In particular, we have al¬ 
ready mentioned that the Weibel instability plays a crit¬ 
ical role in the development of the tearing mode. There¬ 
fore, it is expected that the mirror instability can also 
be important in the stability properties of the tearing 
mode. Indeed, it has already been identified that the 
force responsible for both anisotropic tearing and mirror 
instability have very similar origins^^’®®. And since we 
saw a correlation between stabilization of tearing mode 
and formation of bifurcated structure in the CIC runs, it 
is worthwhile to analyze the relation between bifurcation 
and mirror instability happening as a consequence of an 
imposed temperature anisotropy. In order to do that, let 
us investigate the time evolution of the quantities related 
with bifurcation at different distances from the center of 
the current sheet. The results are shown in Fig. 19. 






t Qei t Qei 

Figure 19. (Short) time evolution for some quantities related 
with bifurcation, averaged along y, at different distances x 
from the center of the CS, for the run TSC-360ppc-A1.96. a) 
Electron density rie- b) Electron bulk speed 14,2. c) Electron 
current density Je,z- d) Ion current density Ji, 2 . We do not 
show the ion density rii since its values in the course of the 
evolution are always very close to Ue, i.e.: the quasineutrality 
condition is kept everywhere to a large extent. 


We can see that the region close to the center of the 
CS, unstable to the mirror instability from tOci ^ 0.6, 
also develops a drop in both density Ue and the out-of- 
plane electron bulk velocity 14 , 2 . This leads to a de¬ 
crease in the out of plane current density J^. The drop 
in density due to an increase in the magnetic pressure 
is a known signature of the mirror mode®®, that can be 
explained in a fluid approach (the fluctuations in no and 
Bq should be in anti-phase). But the drop in 14,2 is due 
to a different process: the pitch angle scattering that is 
active during that time. The bifurcation is stopped when 
the anisotropy is reduced, and so is the scattering. This 
can be proven by means of a comparison of the previous 
Fig. 18 with Fig. 5. The latter shows that bifurcation 
obtained by the CIC scheme increases as the time goes 
by due to the corresponding increase of the electron tem¬ 
perature anisotropy. But in that CIC case, the scattering 
is not only of pitch-angle type, but also due to the overall 
heating with the associated increment in the total parti¬ 
cle energy. 

It is of central importance to note that even after the 
stabilization of the mirror mode, both Ue and I 4 ,z remain 
with the minimum values reached during the brief time 
of strong mirror interactions, even for very long times. 
The mirror structures are not static, but they oscillates 
around the center of the CS. The consequence of this 
process is the formation, after the mirror structures dis¬ 
appears, of a narrow strip region where the current den¬ 
sity Jz is reduced. The neighborhood of this region that 
never became mirror unstable preserves the original value 
of Jz- As as final result, we observe a bifurcated cur¬ 
rent structure: a drop of Jz at the center of the CS (see 
Fig. 18). We can even explain to what extent the bifur¬ 
cation will modify the structure of the current density. 
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since its maximum extension has to be proportional to 
the component of the mirror wavevector in the parallel 
direction fey. Assuming a fixed fc_L (see explanation in 
Sec. IV C), this quantity has to be proportional to a/A 
by Eq. (23). Then, using Eq. (22) and (18), it turns out 
that this the extension across x of the mirror structures 
is directly proportional to the instantaneous value of the 
anisotropy A~^ — 1, and so is the extension of bifurcation 
(as can be expected). Inside of the mirror structures, the 
average field aligned temperature anisotropy never devi¬ 
ates too much from the marginal stability condition, and 
this explain the bifurcation happening only very close to 
the center of the CS. 


F. Influence of the number of macroparticles per cell on the 
initial Weibel/mirror instabilities 


In Fig. 11, we can see a difference in the develop¬ 
ment of the global anisotropy by using different number 
of macroparticles per cell, especially in the long term 
evolution. The initial stage seems to be similar, but 
it has a critical difference. This can be shown, for the 
case with lowest number of macroparticles and highest 
anisotropy TSC-40ppc-A1.96, through plots similar to 
those of Fig. 12. The results are in Fig. 20. 






Figure 20. Time history of the out of plane magnetic field 
and temperature anisotropy averaged along y, at different dis¬ 
tances X from the center of the CS, for the run TSC-40ppc- 
A1.96. a) and b) Magnetic field \Bz\. c) and d) Electron 
temperature anisotropy The left plots a) and c) 

are for the early stages of the system while the right plots b) 
and d) are for the long term evolution. To be compared with 
Fig. 12. 


We can see that the initial electron temperature 
anisotropy T^^y/T^^x is relaxed from the very beginning: 
there is no the ’’plateau” between 0 < tVld < 0.5 as 
seen in the Fig. 12(c) (for the quantities at the center 
of the CS a; = 0). As a result, the free energy stored 
in the anisotropy can be converted from the very begin¬ 
ning in magnetic energy due to the Weibel instability. 


at a faster rate than with 360ppc (compare Fig. 12(a) 
with Fig. 20(a)). For the points close to the center of the 
CS, this can be attributed to the enhanced noise level 
for 40ppc. This noise can provide a higher initial seed 
magnetic field from which the Weibel generated magnetic 
fields can grow faster. In this context, it is interesting to 
note that the asymptotic value reached by Bz for very 
long times is independent of the distance to the cen¬ 
ter of the CS, being only dependent on the number of 
macroparticles per cell. 


The consequence of the enhanced growth of magnetic 
field Bz in the early stages (for this case TSC-40ppc- 
A1.96) is a faster generation of the conditions for the 
triggering of the mirror instability at the center of the 
CS. We checked that the temperature anisotropy aligned 
with this generated magnetic field reaches Tg ||/Te^_L < 1 
earlier than for the case with TSC-360ppc-A1.96. How¬ 
ever, the Weibel and mirror instabilities saturate more 
or less at the same time, once this initial anisotropy is 
exhausted. This implies that the center of the CS expe¬ 
riences more time unstable to the mirror instability, and 
that can explain why bifurcation is stronger (deeper drop 
in Jz) in the cases with less number of macroparticles per 
cell. This can be shown by comparing the time evolu¬ 
tion of bifurcation related quantities TSC-360ppc-A1.96 
in Fig. 19 with the corresponding plots for TSC-40ppc- 
A1.96 shown in Fig. 21. 
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Figure 21. (Short) time evolution for some quantities related 
with bifurcation, averaged along y, at different distances x 
from the center of the CS, for the run TSC-40ppc-A1.96. a) 
Electron density rie. b) Electron bulk speed 14,2. c) Elec¬ 
tron current density Je,z- d) Ion current density Ji,z. We do 
not show the ion density rii since its values in the course of 
the evolution are always very close to rie, i.e.: the quasineu¬ 
trality condition is kept everywhere to a large extent. To be 
compared with Fig. 19 
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G. Influence of ions in the temperature anisotropy driven 
instabilities 

In order to prove that stabilization of tearing mode 
and bifurcation are only due to electron temperature 
anisotropy, we also tested the effects of ion temperature 
anisotropy by means of two other sets of anisotropic sim¬ 
ulation runs (results not shown here). In the first set, we 
isolated the effects of electron temperature anisotropy by 
imposing initially isotropic ions Ai = 1 and anisotropic 
electrons ^ \ (with A^, € [0.64,1.96]), while all 
the other parameters were identical to the isotropic run 
TSC-40ppc. We found practically the same behaviour as 
shown in Fig. 10 with both anisotropic species: the exis¬ 
tence of a threshold for the stabilization of the Weibel in¬ 
stability and the subsequent generation of mirror modes 
structures and bifurcation (associated with the stabiliza¬ 
tion of tearing mode). The values are very similar to 
those ones found previously, but not exactly the same. 
This is because our choice of parameters implies an ini¬ 
tial ratio Ti ||/Tej| ^ 1, that favours different processes 
related with the Landau damping: the damping of some 
wave modes depends on n/Teji, as well as the thresh¬ 
olds and growth rates of Weibel/mirror instabilities. 

For the second set of anisotropic runs, we isolated 
the effects of ion temperature anisotropy by imposing 
Ai ^ 1 and isotropic electrons Ag = 1 (with Ai S 
[0.64,1.96]), while all the other parameters are identical 
to the isotropic run TSC-40ppc. As expected, we found 
that the ion temperature anisotropy does not change 
too much in the course of a typical simulation, since 
their time scales are much longer than electron time 
scales. The CS only develops a slow growth of tearing 
mode for the same timescales as the isotropic run TSC- 
40ppc. There is only a slow trend towards isotropization 
of temperatures, that can be entirely explained by the 
pitch angle scattering caused by an insufficient number 
of macroparticles per cell. As can be expected, the CS 
success to bifurcate only for values of Ai much larger 
than those found for Af.. Regarding the tearing mode, 
this fact agrees with the theoretical prediction^®: an ion 
temperature anisotropy will stabilize the tearing mode 
growth only if it is very large. This is, however, unlikely 
to be caused by numerical heating, which mostly affects 
the electrons. 

The absence of anisotropy driven instabilities in the 
case with only anisotropic ions is because the Weibel 
instability requires mostly an electron anisotropy to be 
triggered, rather than ion anisotropy. Therefore, with 
isotropic electrons there is no generated magnetic field 
Bz that can provide the necessary conditions for the trig¬ 
gering of the the mirror instability: a field aligned tem¬ 
perature anisotropy due to magnetization and rotation 
of the magnetic field at the center of the CS. As a conse¬ 
quence, there is no bifurcation of the CS at all. Note that 
if we would have magnetic fields as those provided with 
Weibel instability, but with only anisotropic ions, a bifur¬ 
cated structure should also be produced. This is because 


the mirror threshold also depends on the ion anisotropy, 
in more or less the same extent compared to the electron 
anisotropy (see Eq. (18)). 

H. Enhancement of tearing mode with the opposite 
anisotropy 

Now, let us discuss the run TSC-40ppc-A0.64 with an 
initial electron temperature anisotropy Ag = Teji/Te.u = 
0.8 < 1 in the opposite direction to the previously dis¬ 
cussed. In this case, the tearing instability is faster 
than for the isotropic temperature case. The explosive 
phase of the tearing instability takes place already at 
t ^ 20 . This explains the sudden growth of the tem¬ 

perature and anisotropy seen in Fig. 10 for this run: the 
quantities Tg H, and Tg || go off the scale beyond 

t > 2012“^. From the faster merging of magnetic islands 
follows that the tearing mode is enhanced for tempera¬ 
ture anisotropies Ag < 1, in agreement with Ref. 13. 

All the results shown in this section confirm that the 
electron temperature anisotropy Ag can lower the growth 
rates of the tearing mode instability or even stabilizes it 
completely if it is strong enough. And it may also cause 
bifurcation due to the always non negligible noise and as¬ 
sociated numerical scattering present in PIC simulations 
of this physical system. 

V. INFLUENCE OF NUMERICAL PARAMETERS ON THE 
TEARING MODE INSTABILITY 

A. Influence of the shape function 

Tearing mode magnetic islands, regions with closed 
magnetic field lines, are formed very early in the CS evo¬ 
lution. They are due to the initial numerical noise and 
are ordered in a long chain along the CS in y direction. 
Afterwards, those structures start to merge via a coales¬ 
cence process due to the attraction of the current pinches, 
generating larger islands (up to widths of the order of the 
ion inertial length di = cjujpi, marginally stable accord¬ 
ing to the linear theory). This is a characteristic signa¬ 
ture of the tearing instability (see Fig. 3(b)). The evolu¬ 
tion of this process can be estimated by integrating the 
Fourier modes of the vector potential \Az{x,ky)\ in the 
electron singular layer thickness ±/S.ms = V'^PeL In 
this region, most of the electrons are unmagnetized and 
experience meandering instead of gyromotion®. With 
our choice of parameters, seven linearly unstable tear¬ 
ing modes satisfying kyL < 1 are possible (shorter wave¬ 
length modes are evanescent). The time history of those 
modes for the TSC-40ppc run are shown in Fig. 22. In 
that figure we also overplotted the sum of all Fourier 
modes (gray) and the sum of the first seven ones (black). 
Obviously, there is a discrepancy between both curves at 
early times. This means that modes with shorter wave¬ 
lengths than mode number seven {kyL > 1) contribute 
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in a significant way to the total power only in the ini¬ 
tial transient stage. Later on, they can be completely 
neglected. 



Figure 22. Time history of the first seven Fourier modes 
\Az{x,ky)\dx for the run TSC-40ppc. The dashed 
gray line is obtained as a sum of all Fourier modes. The 
dashed-dotted black line is the sum of the first seven Fourier 
modes shown in this plot. 



Figure 23. Time history of the first seven Fourier modes 
\Az{x,ky)\dx for the run CIC-40ppc. The dashed 
gray line is obtained as a snm of all Fourier modes. The 
dashed-dotted black line is the sum of the first seven Fourier 
modes shown in this plot. 

As expected, all tearing modes in Fig. 22 reach the sat¬ 
urated stage at the same time as the reconnected flux (cf. 
Fig. 4(b)). But the time history of the Fourier modes for 
the run CIC-40ppc, Fig. 23, reveals a complete suppres¬ 
sion of the tearing modes, correlated with the respective 
non-growing reconnected flux (see Fig. 4(a)). This is due 
to the stabilizing effect of bifurcation^^’*® that dominates 
the whole structure of the CS. It is also interesting to no¬ 
tice that short wavelength modes with kyL > 1 carry a 
signihcant amount of power for most of the evolution of 
the CS, in contradiction with the linear theory. From 
this, we conclude that it is essential to compute growth 
rates of the tearing instability with the TSC shape func¬ 


tion instead of the often used CIC scheme, unless a very 
large number of macroparticles per cell is used. 


B. Comparison with linear theory 


In order to compare the linear growth rates of the 
tearing mode with the linear predictions, it is necessary 
first to identify the linear stage of the instabilities in the 
simulation results. This is not straightforward in multi- 
mode tearing, because a complex transfer of energy be¬ 
tween different modes takes place at different times. We 
chose the following criterion: the linear regression is per¬ 
formed up to the time at which it starts to be noticeable 
a fast growth of the short-wavelength modes kyL > 1 
(M > 7) that should be stable according to the linear 
theory (we did not show those modes with practically 
vanishing growth rates in Fig. 22). For later times, the 
system enters into a non-linear stage, and thus a com¬ 
parison with linear theory is not possible anymore. A 
comparison between linear and non-linear stages can be 
seen in Figs. 24 and 25. In those figures, the time his¬ 
tory of the vector potential along the center of the CS, 
Az{x = 0), is displayed up to the limit of the linear stage 
{t ^ 35 0“^) and up to the start of the explosive phase 
of reconnection (t ~ 6011“^), respectively. 
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Figure 24. Time history of the profiles of vector potential 
Az{x = 0,t) up to the limit of linear stage of tearing mode 
evolution (tflci 35), for the run TSC-40ppc. 
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Figure 25. Time history of the prohles of vector potential 
Az{x = 0 ) up to the explosive phase of reconnection {tQci ~ 
60), for the run TSC-40ppc. 

Fig. 24 shows a clear dominance of four magnetic is¬ 
lands (their O points are the maxima of and the 
X points the minima) at the end of the linear stage. 
This is in agreement with the linear theory*®: the most 
unstable mode should have M = 4 or kyL ~ 0.545 
(green line in Fig. 22), implying that the dominant is¬ 
lands should have a size Hg in y direction of about 
Hs/{c/LOpi) = 27 r/(fcy c/LOpi) = 6.66 (see, e.g., chp. 7 
in Ref. 90, chp. 12 in Ref. 91, Ref. 31). However, the 
reconnected flux between the X and O points of each one 
of those magnetic islands is very low: only marginally ex¬ 
ceeds the numerical noise level. So, magnetic reconnec¬ 
tion it is not energetically important during this stage. 

Later, in the nonlinear stage, a transfer of power to 
the long wavelengths modes takes place. As can be seen 
in Fig. 25, at the end the whole system is dominated by 
only one tearing island, with a much larger reconnected 
flux than during the linear stage. This is just before 
of the explosive phase of reconnection, when the entire 
structure of the CS is disrupted. 



kyL 


Figure 26. Dots connected by the solid line: simulated growth 
rates 7 of the Fourier modes \Az{x, ky)\ dx (shown in 

Fig. 22) vs wave number kyL for the run TSC-40ppc (see text 
for the calculation method). Dashed line: analytical estimate 
of the growth rate according to Eq. (24) for thin CS and 
rrii = me. 


For the previously described linear stage, we compute 
the growth rates 7 of the tearing instability as function 
of kyL by means of a linear regression of the integrated 
Fourier modes \^z{x,ky)\dx shown in Fig. 22. 

The results are shown in Fig. 26. Some modes with 
kyL > 1 are also included, in order to prove that those are 
linearly damped. Although the criterion for choosing the 
linear stage provides a first order approximation, is not 
good for all modes since some of them reach a non-linear 
stage before. For that reason, we chose several windows 
for performing the linear fit from very early and up to 
t ~ 35fl“/. Each one provides a different growth rate, 
being closer to each other for a exponential growth, and 
we plot the average of those values in Fig. 26. The error- 
bars are the standard deviation of those values. In the 
Fig. 26, the dashed lines are from an analytical formula 
for thin CS {pi/L ~ 1, Ref. 89): 

7 2'/2 / kyL(2 + kyL){l — kyL) 

^ci 'JZ \l) 1 + a (^)^ 

=0.2591 kyL (2 -|- kyL){l — kyL). 

This expression for the collisionless ion tearing mode 
growth rate was derived within the Vlasov formalism 
by matching internal and external solutions for the per¬ 
turbed vector potential with respect to the electron sin¬ 
gular layer thickness A.^s- Note that this approxima¬ 
tion to the tearing mode growth rate considers equal 
masses for both ion and electrons mi/mg = 1. Nu¬ 
merical solutions of the linearized Vlasov equation^’®^’®* 
have shown that for realistic ion to electron mass ra¬ 
tios {rrii/me = 1836), the growth rates are reduced by 
a factor of 1.5 to 2 (see section 10.4 of Ref. 94). For 
the intermediate mass ratios used in our simulations, 
as expected, the growth rates are reduced by a smaller 
amount. This explains, in part, the difference between 
the maximum value of the growth rate predicted by the 
analytical formula 'y/ilci = 0.163 for the most unstable 
mode kyL ~ 0.545. But the main reason for the discrep¬ 
ancy in the overall simulated curve with the theoretical 
one is due to the coupling between the unstable modes 
for the multimode tearing, as it has been reported by 
previous works (see, e.g.: Ref. 42). 

C. Influence of the number of macroparticles on tearing 
mode growth rates 

In spite of all the previous comparison with the lin¬ 
ear theory, it is important to remark that the dominant 
modes with higher power, as well as the values of the 
growth rates, vary with the number of macroparticles 
used in the simulation. There are two main reasons for 
this: many non-linear effects absent in the linear theory 
are enhanced by the numerical noise. Also, the domi¬ 
nance of the fastest growing Fourier modes is delayed by 
the same reason®®. Thus, it is expected to have a better 
match with the theoretical value for a larger number of 




















21 


macroparticles, i.e.: lower noise level. However, in agree¬ 
ment with the linear theory, the most unstable mode with 
higher growth rates, M = 4, is always the same. In or¬ 
der to prove this claim, we run an additional simulation 
denoted as TSC-640ppc, with the same parameters as 
TSC-40ppc but with 16 times more macroparticles per 
cell. The latter reduces the initial noise to 1/4 of its 
original value. The results for the time history of the in¬ 
tegrated Fourier modes are shown in Figs. 27, while the 
respective growth rates are in Fig. 28 (calculated with 
the same previously explained method). 



Figure 27. Time history of the first seven Fourier modes 
\Az{x, ky) \ dx for the run TSC-640ppc. To be com¬ 
pared with Fig. 22, for the run TSC-40ppc with sixteen times 
less macroparticles per cell. The dashed gray line is obtained 
as a sum of all Fourier modes. The dashed-dotted black line 
is the sum of the first seven Fourier modes shown in this plot. 



kyL 


Figure 28. Dots connected by the solid line: Simulated growth 
rates 7 of the Fourier modes \Az{x, ky)\ dx (shown in 

Fig. 27) vs wave number kyL for the run TSC-640ppc (see text 
for the calculation method). Dashed line: analytical estimate 
of the growth rate according to Eq. (24) for thin CS and 
rtii = rrie. Compare with the run TSC-40ppc (sixteen times 
less macroparticles per cell) shown in Fig. 26. 

Comparing the evolution of the Fourier modes in 
Fig. 27 and Fig. 22, one sees a clear dominance (higher 
power) of the mode M = 4 {kyL ^ 0.545) in case of more 


macroparticles per cell. A comparison for the growth 
rates in Fig. 28 and Fig. 26 shows that this same mode 
is the most unstable in both cases. However, the growth 
rates for TSC-640ppc are smaller due to the lower numer¬ 
ical noise. Note also that the growth rates for the modes 
M=1 and M=7 have much smaller error bars than for 
the run TSC-40ppc (compare with Fig. 26). This implies 
that the reduction of the numerical noise affects mainly 
those modes by extending the duration of their linear 
growth phase. 

VI. DISCUSSION AND CONCLUSIONS 

Before drawing the conclusions from this work, it is 
necessary to discuss two results regarding diagnostics and 
simulations which are not shown here. 


A. Influence of the boundary conditions 

We used conducting boundary conditions in the x di¬ 
rection in our simulations. It has been claimed that those 
boundaries are not a realistic choice for applications to 
space plasmas. To clarify this issue, we investigated their 
influence on our simulation results. First, we detected 
waves generated at the center of the current sheet in the 
beginning of the simulation caused by the initial numer¬ 
ical noise. These waves propagate outwards becoming 
reflected back and forth between the conducting bound¬ 
aries and the center of the current sheet. This leads leads 
to a periodic exchange of magnetic and kinetic energy. 
The generation of such magnetosonic waves was already 
reported before (Refs. 96 and 97). They, however, do not 
essentially affect the physical results. We verified this 
conclusion by changing the simulation box length across 
the reflective direction x (results not shown here). 

B. Influence of the mass ratio 

We found that the numerical anisotropic heating is es¬ 
pecially critical for higher mass ratios, such as mi/mg = 
180 and above, compared to the often used reduced mass 
ratio 25. The difference is due to the separation of the 
time scales of ions to electrons: Indeed, while 

time scales of the evolution of the tearing mode instabil¬ 
ity depends only on the number of iterations per¬ 
formed by the code is proportional to uip^. This has 
two consequences. First, the number of iterations before 
the onset of the same physical phenomenon is smaller 
for lower mass ratios (since oc And 

second, since stochastic heating increases the kinetic en¬ 
ergy linearly with the number of timesteps (see section 
9.2 of Ref. 20), the accumulation of numerical errors is 
less important for lower mi/me. 

In order to prove this claim, we ran simulations (not 
shown here) with the same parameters but lower mass 
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ratio (mi/me=25), and, thus, a smaller ratio 
We found no significant numerical heating with a CIC 
scheme. As expected, there is no important differ¬ 
ence between the evolution of energy and temperatures 
anisotropies with either TSC or CIC shape function. This 
is the reason why no bifurcation was seen in the past PIC- 
code simulations, when lower mass ratios were used. In 
addition to that, growth rates of the tearing mode do 
not show important differences between runs using dif¬ 
ferent schemes. Nevertheless, it is known that the choice 
of too small mass ratios can reveals misleading results^^. 
Therefore, the choice of higher order shape functions be¬ 
comes important when more realistic results (for higher 
mass ratios) are desired. 

The mass ratio is also a key parameter for the ini¬ 
tial evolntion in the initially anisotropic rnns discnssed 
in Sec. IV. Indeed, both Weibel and mirror instabilities 
have growth rates that depends on the mass ratio nrii Inie 
(at least, in the numerical solutions for the fully kinetic 
expressions not shown here). By choosing more realistic 
valnes, the growth rates are in general increased. This 
implies that a PIC simulation with a too low mass ratio 
will underestimate the importance of those instabilities, 
leading to an incorrect physical evolution of the system. 
Moreover, the dependence of those kinetic instabilities 
(as well as many other ones) on rm/me is in general dif¬ 
ferent. Therefore, the change of the mass ratio would 
may also imply a change in the dominant instability of 
the system, affecting dramatically the physical predic¬ 
tions of a PIC simulation (see Ref. 58 and 98 for further 
details). 


C. Influence of the restriction to an antiparallel magnetic 
field configuration 

It is important to remark that all the effects related 
with the anisotropy mentioned in this paper are obtained 
for an antiparallel magnetic field configuration, neglect¬ 
ing any finite (guide) magnetic field in the out-of-plane 
direction. Under the influence of a (strong enough) guide 
field, electrons become magnetized (gyrotropic) in this 
gnide field with a highly isotropic temperature in the re¬ 
connection plane. As a result, the stabilization and satu¬ 
ration of the tearing mode in the guide field case will not 
depend on the electron temperature anisotropy. Instead 
of that, it has more to do with processes like the electron 
trapping in magnetic islands, etc. (see, e.g.. Refs. 31 and 
42). 

Therefore, the proper modelling of the electron tem- 
peratnre anisotropy is especially critical for the antipar¬ 
allel field case, while in case of a finite guide field, other 
quantities are more important for the stabilization of the 
tearing mode, such as the resolution of the electron Lar- 
mor radius on the guide field. 


D. Conclusions 

By means of PIC simulations of thin collisionless Har¬ 
ris current sheets, we have shown that higher order 
macroparticle shape functions allow to obtain more real¬ 
istic (less numerically determined) PIC code results. In 
particular, the use of the TSC scheme results in a much 
better energy conservation over long times than the tra¬ 
ditionally nsed CIC shape function. In an efficient way, 
the TSC shape function suppresses numerical collisions 
and anisotropic heating better than the CIC shape func¬ 
tion, i.e., at less compntational cost than a larger nnmber 
of macroparticles per cell. This is in agreement with the 
results of a previous study of laser wakefield interactions 
by PIC code simulations^^. We fonnd that the choice of 
a TSC shape function is especially important for larger 
(more realistic) mass ratios and for long PIC rnns. 

Fnrthermore, lower order shape functions (such as the 
often used CIC) or a too small number of macroparticles 
per cell enhance the numerical collisions inherent to all 
PIC codes, leading to non-physical results. Since numer¬ 
ical collisions correspond to irreversible processes, they 
enhance the entropy of the system. We found that for 
higher order shape functions (TSC instead of CIC), the 
increase of entropy slows down more than by the use of 
a larger number of macroparticles per cell (see Fig. 7). 
This is because the diffusion coefficient of the effective 
Boltzmann equation for the macroparticles, as a mea¬ 
sure of the collisionality, depends explicitly on the shape 
function (see Eqs. (B2) and (Bl)). 

Numerical collisions correspond to an effective scatter¬ 
ing. They cause unphysical numerical heating^^, leading 
to an artificial increase in the total energy of the system 
(see Fig. 2(a)). Note that the heating due to effective 
scattering is different from grid heating^^. The latter de¬ 
pends mainly on the number of macroparticles per cell. It 
becomes dominant when the grid size is larger than the 
Debye length, i.e.: not in our case. Numerical heating 
affects mostly electrons. This in agreement with previ- 
ons simulations and analytical estimations (See section 
9.2 of Ref. 20). It was shown that the numerical heating 
is inversely proportional to the mass of the macroparti¬ 
cles considered, being thus negligible for the ions. It is 
also important to remark that the heating time depends 
strongly on the type of shape function, while the collision 
time does not (it only has a weakly dependence on the 
shape function. See Eq. 9.22 of Ref. 20 or Ref. 33 for 
the 2D electrostatic case). Following Ref. 33, the heating 
time measures the energy losses in an ideal homogeneous 
plasma without external electromagnetic fields, ideally 
being infinite in that system. In our case, we expect 
physical heating of the plasma only at later (non-linear) 
stages of the tearing mode evolution, due to particle ac¬ 
celeration caused by the reconnected electric field in the 
ont-of-plane direction. 

The observed electron heating is anisotropic, as one 
can see in Fig. 6. This was already noticed in an early 
work®^ and has to do with the anisotropic nature of col- 
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lisions in 2D3V PIC code simulations in magnetic fields. 

It is known that temperature anisotropies can affect 
the stability of a current sheet. Indeed, we confirmed the 
theoretical anisotropy threshold for the stabilization of 
the tearing mode instability^^, given by Eqs. (9) or (8), 
by imposing anisotropy and by letting it develop in simu¬ 
lations with low order shape functions or too small num¬ 
ber of macroparticles per cell. 

We also confirmed that the temperature anisotropy 
drives additional micro-instabilities. For strong enough 
(electron) temperature anisotropies, as they develop in 
simulations with low order shape functions (CIC), prone 
to numerical heating, current sheets bifurcate by a reduc¬ 
tion of the out-of-plane component of the current density 
Jz around their center, which corresponds to a reduc¬ 
tion of the out-of-plane component of the electron bulk 
velocity 14 ,z- In agreement with previous studies, we 
confirmed that bifurcated current sheets inhibit the de¬ 
velopment of the tearing mode instability. Bifurcation 
of Harris-type current sheets occurs around near their 
center, where the magnetic field vanishes. It can be due 
to the electron chaotic scattering^® after anisotropies are 
enhanced by, e.g., numerical collisions in electromagnetic 
PIC code simulations. Note that these collisions are in¬ 
versely proportional to the magnetic field strength®®. 

We proved that a Harris current sheet can sponta¬ 
neously bifurcate even in schemes with good energy con¬ 
servation (TSC shape function), in case of an initializa¬ 
tion with a high enough level of electron temperature 
anisotropy A^.. Only in the limit of negligible scattering, 
current sheets do not bifurcate. For our parameters, the 
threshold for the onset of bifurcation is around A^. > 1.2, 
just a little bit above the analytical estimate for tearing 
mode stabilization given by Fq. (8) (which is related with 
the Weibel instability criterion). Although this threshold 
was derived for a collisionless Vlasov plasma, it can still 
be useful for predicting bifurcation for levels of numerical 
collisions feasible to find in typical PIC simulations. 

We also clarified the mechanism of the initial tempera¬ 
ture relaxation anisotropy seen in some of the cases with 
an initially imposed electron temperature anisotropy. Al¬ 
though the initial state with a temperature anisotropy Ae 
is still a Harris equilibrium, under the right conditions it 
can be unstable, in addition to the tearing mode, to the 
Weibel instability close to the unmagnetized center of the 
current sheet. A high enough initially imposed electron 
anisotropy (see threshold condition Eq. 14) will necessar¬ 
ily have to be relaxed shortly after the initial stage via 
Weibel instability, at much shorter timescales than the 
ones associated with the slow growth of tearing mode. 
This instability generates filamented structures of the 
out-of-plane magnetic field relaxing the anisotropy in 
addition to magnetize the center of the current sheet. In 
addition, this magnetization provides the conditions for 
the triggering of the mirror instability afterwards. Both 
instabilities grow together until the saturation stage, but 
the regions around the center affected by the mirror in¬ 
stability will develop a bifurcation of Jz- The details of 


the evolution depends on the number of macroparticles 
per cell, since the growth of the Weibel generated mag¬ 
netic field has to start from the magnetic field originated 
from the initial numerical noise, but this physical process 
should always be observed under the right conditions. 

By isolating the effects of initially imposed electron or 
ion temperature anisotropies, we confirmed that only the 
electron temperature anisotropy Ag is important for the 
stabilization of the tearing mode, even if there is an ion 
temperature anisotropy Ai 0: the latter does not play 
a role in this process, at least for the small values of A^ 
analyzed here. An ion temperature anisotropy can stabi¬ 
lize a tearing mode only if Ag is practically zero and Ai 
is very large^®. Moreover, the fast relaxation of an ini¬ 
tially imposed temperature anisotropy via Weibel/mirror 
instabilities does not take place with only Ai ^ 0. This is 
because the Weibel generated magnetic field that provide 
the conditions for the mirror modes, and the associated 
bifurcation later on, requires an electron anisotropy to be 
triggered. 

Finally, taking into account all those considerations, 
we found growth rates of the simulated tearing instabil¬ 
ity that match the theoretical predictions (see Figs. 26 
and 28). There are discrepancies, however, that are pro¬ 
duced by the nonlinear exchange of energy among dif¬ 
ferent modes allowed in a multimode tearing situation, 
and to a lesser degree by the enhanced noise level (coarse 
graining effect), due to the reduced number of physical 
particles represented by one macroparticle in a PIC code. 
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Appendix A: Shape function 

The use of a spatial grid in PIC codes implies that ev¬ 
ery macroparticle should be distributed in the real space, 
as described by a shape function S'(x)^®’^°. Shape func¬ 
tions distribute the charge of the macroparticles over the 
grid and determinate the action of the electromagnetic 
forces on the macroparticles. S(x) can be defined by the 
following expression for a single-macroparticle distribu¬ 
tion function, 

f{x,v)=R6{v-Vn)S{x-Xn), (AI) 

where Xn is the position around which the macroparticle 
is distributed (smeared out), fT„ is its velocity and R is the 
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ratio of the number of physical particles to (numerical) 
macroparticles. It is appropriate to normalize S{x) to 1 
and choose it separable: S{x — Xn) = S{x — Xn)S{y — 
yn)S(z — Zn)^^ ■ Eq. (Al) implies that the macroparticles 
have a hnite spatial extent. It can be characterized by the 
weight function W (x), which is obtained by integrating 
the shape function over the cell volume 

/•Xc+Aa:/2 

W(Xc—Xn)= / S{x' — Xn)dx' . (A2) 

J Xc —Aa:/2 

Hence, W{x) characterizes the overlap of the spatially 
distributed macroparticle and the grid cell of size Ax 
centered around the mesh point Xc- 

In practice, only a small number of different shape 
functions are used. Higher order functions are generally 
more accurate (since they reduce the aliasing effects as¬ 
sociated with the undersampling by the interpolation®®). 
On the other hand, they are also computationally more 
expensive. The use of higher order shape functions re¬ 
duce the error in the calculation of non-physical forces 
between the macroparticles smeared over the grid. Thus, 
the corresponding PIC simulated plasma behaves more 
similar to a real plasma^®®. 

The simplest shape function is the zero-order weighting 
scheme called NGP (Nearest Grid Point), with the (one¬ 
dimensional) weight function 


W{x) 


1 ife<l/2, 

0 otherwise. 


(A3) 


given by: 

\\-e ife<iA 

ifI/2<e<3/2, (A5) 

I 0 otherwise. 

A TSG weighting scheme is much smoother than the GIG 
scheme. It is, however, computationally more expensive 
and therefore not so often used in PIG codes calcula¬ 
tions. A third order (cubic) scheme is the PQS (Piecewise 
Quadratic Spline), with a weight function: 

r 1 (4 - se) if 0 < e < 1, 

Wix) = A (2 - ^3) if 1 < e < 2, (A 6 ) 

I 0 otherwise. 

These different shape functions and their associated 
weight functions are illustrated in Fig. 29 for ID and 
in Fig. 30 for 2D. 



X/AX X/AX 


Figure 29. Left panel: ID shape functions S{x). Right panel: 
ID weight functions W{x). Continuous line: NGP. Dashed 
line; CIC. Dotted line: TSC. Dash-dotted: PQS. Those fig¬ 
ures are for a macroparticle located at the origin Xn = 0. 


Here, ^ = |x — Xn\/Ax is the relative distance with re¬ 
spect to the center of the macroparticle. A NGP func¬ 
tion weighting assigns all the charge of a macroparticle 
to its nearest grid point. The corresponding S{x) is a 
Dirac delta function: each macroparticle is concentrated 
at one location. Although fast, the NGP scheme is rarely 
used since it is very noisy (due to the enhanced numerical 
scattering among the macroparticles). 

The next higher order weighting scheme can be ob¬ 
tained through a convolution of the NGP weighting func¬ 
tion with itself (the so-called b-splines are obtained by 
applying this method successive times. See, e.g. Ref. 21). 
This CIC (Cloud in Cell) scheme, the most used one in 
PIC codes, is defined by the weight function: 
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W{x) 


1-^ iA<l, 

0 otherwise. 


(A4) 


The CIC scheme distributes the charge of every 
macroparticle between the two nearest grid points by 
means of a linear interpolation. It provides a good com¬ 
promise between the smoothness of the interparticle force 
and computational speed. 

Next, there is a second order scheme called TSC (Tri¬ 
angular Shaped Cloud), with a quadratic weight function 


Figure 30. 2D weight functions W(x)W(y) for different 
weighting schemes, a) NGP (Nearest Grid Point), b) GIG 
(Gloud in Gell). c) TSG (Triangular Shaped Gloud). d) PQS 
(Piecewise Quadratic Spline). 


There are three consequences of choosing successive 
higher order shape functions. First, an increase of the 
width of W{x) by Ax. Second, a smoother continuity^^ 
of the charge density distribution and the electromag¬ 
netic fields. Those functions increase their continuity 
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class from C" to with n the order of the interpola¬ 

tion scheme. The third and most important consequence 
for the goal of this paper is the reduction of numerical 
collisions between the macroparticles by using higher or¬ 
der shape functions®^. See next Appendix B for their role 
in PIC codes. 


Finally, it is important to remark that numerical colli¬ 
sions can be interpreted as a scattering of macroparticles, 
and therefore they can cause a non-physical particle heat¬ 
ing (the second numerical heating mechanism discussed 
in Sec. I). 


Appendix B: Numerical collisions 
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Numerical collisions are principally inherent to all PIC 
codes. In fact, due to the calculations on a spatial grid, a 
finite time step and the representation of many physical 
particles by macroparticles (coarse-graining), PIC codes 
do not model a strictly collisionless plasma, which would 
fully obey the Vlasov-Maxwell equations (see Ref. 97 for 
a extended discussion about this issue). Instead of the 
Vlasov equation, PIC codes practically solve a kinetic 
Boltzmann equation for the distribution function, with 
an effective numerical collision operator absent in the 
Vlasov equation. This collision operator can be estimated 
as being proportional to^®: 


dt 


—— ] oc dk 


sHk) 


\e{k,k-v)\^ 

dif S{k ■ V — kp ■ if,ujg). 


E (Bi) 


For more detailed discussion, see appendix E of Ref. 19. 
Here k and K are the finite difference gradient and Lapla- 
cian operators, respectively (associated with k and k^, re¬ 
spectively). kp = k — p - kg, kg = 2 tt /is the grid 
wave number, cug = 2 TTjAt is the characteristic frequency 
of the time stepping, S{uj,ujg) = ~ i® 

a periodic delta-function comb, e is the plasma dielec¬ 
tric function which does not only depend on the physical 
model being analyzed, but also depends on the time in¬ 
tegration scheme, the conservation properties of the al¬ 
gorithm, the shape function and other details of a PIC 
code. 

As one can see in Eq. (Bl), the collision operator in¬ 
volves directly the Fourier transform of the shape func¬ 
tion S{k) as well as the grid size Ax and the simulation 
time step At. It can be rearranged in the form of a 
Fokker-Plank term with effective diffusion Dij and drag 
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The coefficients Dij and Ai depend on the spread of the 
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sion coefficients decrease as well as the influence of the 
numerical collisions. 
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